1/2 


AD-A171  258  THREE-DIMENSIONAL  ANALVTICAL  MODELING  OF 

DIFFUSION-LIMITED  SOLUTE  TRANSPORT(U)  AIR  FORCE  INST  OF 
TECH  HRIGHT-PATTERSON  AFB  OH  M  N  G0LT2  JUL  96 
UNCLASSIFIED  AFIT/CI/NR-86-1310  F/G  8/8  NL 


MICROCOPY  RtSOiUHOM  US! 

NATIONAL  M(Mj  *  •  * 


OTIC  Fils  COPY  AD- A 171  258 


|  REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIOf.-^—^ 
BEFORE  COMPLETING  FORM 

t.  REPORT  NUMBER 

AFIT/CI/NR  86-  131D 

2.  GOVT  ACCESSION  NO. 

3.  RECIPIENT'S  CATALOG  NUMBER 

«.  TITLE  (and  Submit) 

Three-Dimensional  Analytical  Modeling  Of 
Diffusion-Limited  Solute  Transport 

5.  TYPE  OF  REPORT  6  PERIOD  COVERED 

Tlf9^/DISSERTATI0N 

6.  PERFORMING  ORG.  REPORT  NUMBER 

7.  AUTHORS) 

Mark  Neil  Goltz 

8.  CONTRACT  OR  GRANT  NUMBERfaJ 

9.  PERFORMING  ORGANIZATION  NAME  AND  AODRESS 

AFIT  STUDENT  AT:  Stanford  University 

10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  8  WORK  UNIT  NUMBERS 

II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

AFIT /NR 

WPAFB  OH  45433-6583 

12.  REPORT  DATE 

1986 

13.  NUMBER  OF  PAGES 

172 _ 

U.  MONITORING  AGENCY  NAME  A  ADDRESSf/f  different  from  Controlling  Offic e) 

_ 

IS.  SECURITY  CLASS,  (ol  thh  report) 

UNCLAS 

TT  I'TT  mi  rr— 

16.  DISTRIBUTION  STATEMENT  (ol  this  Report)  ■ - i — |—  |  _  - 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED  I—'  1  1 

r^OELECTE 

AU62  8B86 

i 

17.  DISTRIBUTION  STATEMENT  (ol  Iht  abstract  tnltrtd  In  Block  20,  II  dllltrenl  horn  Report) 

_ J _ J 

B 

_ 

18.  supplementary  notes  _  / VA. -  L _ 

CnMnE^WOLAVElf^ 

APPROVED  FOR  PUBLIC  RELEASE:  IAW  AFR  190-1  W&n  for  Research  and ' 

Professional  Development 
AFIT/NR 

19.  KEY  WOROS  (Continue  on  reverse  aide  If  necessary  and  Identify  by  block  number> 

20.  ABSTRACT  (Continue  on  reverse  side  It  necessary  end  Identity  by  block  number) 

ATTACHED. 

DD  1  JAN  73  1473  EDITION  OF  I  NOV  65  IS  OBSOLETE 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  fttfctn  Dele  Entered) 


THREE-DIMENSIONAL  ANALYTICAL  MODELING  OF 


DIFFUSION-LIMITED  SOLUTE  TRANSPORT 


A  DISSERTATION 

SUBMITTED  TO  THE  DEPARTMENT  OF  CIVIL  ENGINEERING 
AND  THE  COMMITTEE  ON  GRADUATE  STUDIES 
OF  STANFORD  UNIVERSITY 
IN  PARTIAL  FULFILLMENT  OF  THE  REQUIREMENTS 
FOR  THE  DEGREE  OF 
DOCTOR  OF  PHILOSOPHY 


By 

Mark  Neil  Goltz 
July  1986 


I  certify  that  I  have  read  this  thesis  and  that  in  my 
opinion  it  is  fully  adequate,  in  scope  and  auallty,  as 
a  dissertation  for  the  degrae  of  Doctor  of  Philosophy. 

r  (Principal  Adviser) 


1  certify  that  I  have  read  this  thesis  and  that  in  ay 
opinion  it  is  fully  adequate,  in  scope  and  quality,  as 
a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


1  certify  that  1  have  read  this  thesis  and  that  in  ay 
opinion  it  is  fully  adequate,  in  scope  and  quality,  as 
a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


I  certify  that  I  have  read  this  thesis  and  that  in  ay 
opinion  it  is  fully  adequate,  in  scope  and  quality,  as 
a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


I 


Approved  for  the  University  Committee 

on  Graduate  Studies: 


Dean  of  Graduate  Studies 


ACKNOWLEDGMENTS 


Until  recently,  completion  of  this  dissertation  seemed  to  be  an 
unattainable  goal.  Now,  with  the  help  of  many,  that  goal  has  been 
reached.  It  is  with  pleasure  and  gratitude  that  I  thank  those  who 
helped  me  attain  what  I  thought  was  unattainable. 

First  and  foremost,  I  am  indebted  to  Paul  Roberts,  my  principal 
adviser,  who  eased  the  frustrating  and  painful  process  of  developing  a 
research  topic,  expedited  the  conduct  of  my  research,  and  assisted 
immeasurably  in  the  creation  of  this  final  document.  His  ideas,  hard 
work,  and  efforts  on  my  behalf  are  sincerely  acknowledged. 

I  would  also  like  to  thank  the  members  of  my  dissertation  commit¬ 
tee,  David  Freyberg,  George  Parks,  and  Martin  Relnhard,  whose  timely 
guidance  helped  keep  my  research  on  track. 

Most  of  my  waking  hours  over  the  past  four  years  have  been  spent  in 
the  company  of  the  "Terman  moles.”  These  denizens  of  the  Terman  Engi¬ 
neering  Center  basement  have  provided  friendship.  Intellectual  stimula¬ 
tion  and  help,  and  needed  distractions.  My  thanks  to  Helen  Dawson,  Lew 
Semprini,  Craig  Criddle,  Tim  Vogel,  Sandy  Robertson,  Kim  Hayes,  Tom 
"Clarence”  Black,  Susan  Stipp,  George  Redden,  Lambis  Papelis,  Gary 
Curtis,  Tong  Xinggang,  Karen  Gruebel,  Meredith  Durant,  Scott  Summers, 
and  Doug  Kent.  Special  thanks  to  Bill  Ball,  who  Introduced  me  to  the 
mysteries  of  the  laboratory.  I  would  also  like  to  especially  thank  my 
good  friend  and  office  mate,  Avery  Demond,  whose  sense  of  humor  kept  me 
going,  critical  reviews  kept  me  thinking,  and  editing  skills  kept  my 
Brooklynese  out  of  this  work. 

Several  people  outside  Stanford  made  significant  contributions  to 
this  work.  Rien  van  Genuchten  provided  the  computer  code  that  was  used 
as  the  basis  for  some  of  the  codes  developed  during  this  research. 
Discussions  and  correspondence  with  A1  Valocchi  proved  to  be  of  great 
assistance.  In  particular,  Dr.  Valocchi 's  help  in  verifying  the  analyt¬ 
ical  solution  presented  herein  was  extremely  valuable.  Doug  Mackay 
provided  much  encouragement  and  insight,  especially  during  the  trying 
early  periods  of  my  research. 


I  would  like  to  express  ay  thanks  to  the  "cast  of  thousands"  who 
were  involved  with  the  Borden  field  experiment.  The  interesting  results 
of  that  experiment  provided  the  motivation  for  this  research.  Thanks  to 
the  Environmental  Protection  Agency,  who  funded  the  Borden  experiment 
and  provided  logistical  support  for  my  research  under  Contract  EPA-CR- 
808851.  Much  thanks  to  the  U.S.  Air  Force,  who  sent  me  to  Stanford  for 
the  best  (and  busiest)  assignment  I  ever  had.  Also,  sincere  thanks  to 
Ditter  Peschcke-Koedt,  who  magically  transformed  my  rough  first  draft 
into  a  clean,  professional  document. 

Lastly,  I'd  like  to  thank  my  family.  Mi  Suk,  my  wife,  bore  the 
burden  of  caring  for  our  boys,  Hugh  and  Eric,  while  I  worked  on  this 
dissertation.  Mi  Suk's  love,  support,  and  patience  were  instrumental  in 
the  successful  completion  of  this  work. 


ABSTRACT 


Previous  experimental  work  conducted  using  laboratory  soil  columns 
has  shown  that  diffusion  into  regions  of  immobile  water  can  have  a  large 
effect  on  solute  transport  through  porous  media*  This  study  focuses  on 
the  development,  analysis,  and  application  of  an  analytical  model  which 
Incorporates  the  diffusion  mechanism  into  the  traditional  three- 
dimensional  advective/dispersive  solute  transport  equation. 

By  consecutively  applying  the  Laplace  transform  in  time  and  the 
Fourier  transform  in  space,  analytical  solutions  are  derived  for  the 
coupled  partial  differential  equations  which  describe  three-dimensional 
advective/dispersive  transport  through  regions  of  mobile  water  and 
Fickian  diffusion  through  immobile  water  regions  of  simple  geometry 
(spherical,  cylindrical,  and  layered). 

To  assist  in  the  analysis  of  the  models,  a  mod  . led  form  of  Aris’ 
method  of  moments  is  presented,  which  permits  the  calculation  of  the 
spatial  and  temporal  moments  of  the  three-dimensional  diffusion  models, 
without  having  to  invert  the  Laplace  or  Fourier  transformed  solutions. 
Using  this  method,  the  moments  of  the  diffusion  models  are  compared  with 
one  another,  with  the  moments  of  a  model  that  assumes  equilibrium  advec¬ 
tive/dispersive  transport,  and  with  the  moments  of  a  model  that  assumes 
a  first-order  rate  law  governs  mass  transfer  between  the  mobile  and 
immobile  regions.  .  The  method  of  moments  is  also  used  to  analyze  the 
differences  in  the aoatial  and  temporal  moment  behavior  of  each  trans¬ 
port  model  under  discuaslon. 

Finally,  the  results  of  a  field  experiment  conducted  to  study 
sorbing  solute  transport  are  presented  and  interpreted  using  these 


models.  It  is  shown  that  the  first-order  rate  snd  diffusion  models 
offer  one  plausible  explanation  of  experimental  observations  which  are 
unexplainable  using  the  traditional  advective/dispersive  model  approach. 
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CHAPTER  1 
INTRODUCTION 


BACKGROUND 

In  recent  years,  the  problem  of  groundwater  contamination  has  re¬ 
ceived  widespread  public  attention.  Hazardous  chemical  waste  is  being 
generated  at  the  rate  of  60  million  tons  annually  (U.S.  EPA,  1980). 
Whether  intentionally  disposed  of  or  accidentally  spilled,  some  of  this 
waste  can  eventually  reach  the  groundwater  and  contaminate  it.  Ground- 
water  may  transport  the  contaminants  from  the  Initial  disposal  site  to 
an  area  where  a  threat  to  public  health  may  be  posed.  Therefore,  it  is 
Important  to  understand  the  processes  affecting  the  transport  of  these 
contaminants  in  the  subsurface  environment.  In  a  study  of  the  problem, 
the  Panel  on  Groundwater  Contamination  of  the  National  Research  Council 
remarked , 

Reliable  and  quantitative  prediction  of  contaminant 
movement  can  be  made  only  if  we  understand  the  pro¬ 
cesses  controlling  transport,  hydrodynamic  disper¬ 
sion,  and  chemical,  physical  and  biological  reactions 
that  affect  soluble  concentrations  in  the  ground. 

(NRC,  1984). 

Traditionally,  the  mathematical  models  used  to  describe  solute 
transport  in  groundwater  flow  systems  have  been  premised  on  the  advec- 
t ion/dispersion  equation.  Advection  refers  to  the  average  motion  of 
solute  due  to  the  groundwater  flow.  Dispersion  describes  spreading  of 
solute  about  the  mean  displacement  position.  Dispersion  is  usually 
attributed  to  molecular  diffusion,  so-called  mechanical  dispersion,  and 
spatial  variability.  Mechanical  dispersion  is  caused  by  local  velocity 
variations  along  tortuous  flow  paths  and  the  velocity  distribution 
within  each  pore  (Bear,  1979).  Spatial  variability  produces  velocity 
variations  on  a  macroscopic  scale.  In  the  advection/disperslon  model, 
the  mechanical  dispersion  mechanism  and  spatial  variability  effects  are 
often  assumed  to  be  diffusive  (that  Is,  the  dispersive  flux  due  to 
mechanical  dispersion  and  spatial  variability  can  be  expressed  by  a 
Fickian  type  law).  Usually,  the  effect  of  molecular  diffusion  is  as¬ 
sumed  negligible  in  comparison  with  the  effect  of  mechanical  dispersion 
(Gillham  et  al.,  1984). 


If  the  solute  sorbs  onto  the  aquifer  material,  it  is  often  assumed 
that  sorption  is  instantaneous,  linear,  and  reversible.  These  assump¬ 
tions  permit  modeling  of  sorbing  solute  transport  using  the  advection/ 
dispersion  equation  as  well  (Bear,  1979). 

The  advection/disperslon  equation  is  usually  implemented  under  the 
assumption  that  the  dispersion  coefficient  is  a  constant  property  of  the 
porous  medium  and  the  mean  velocity.  However,  results  of  laboratory  and 
field  Investigations  have  shown  that  the  dispersion  coefficient  depends 
on  the  scale  of  the  test  or  the  size  of  the  domain  through  which  the 
solute  travels,  thereby  invalidating  the  assumption  underlying  the 
advectlon/dispersion  model  in  its  simple  form  (Gillham  et  al.,  1984). 
Stochastic  models  have  been  used  to  account  for  this  scale  effect.  Sto¬ 
chastic  models  assume  that  the  statistical  structure  of  a  conductivity 
field  in  a  heterogeneous  medium  can  be  estimated.  The  models  also 
assume  that  the  heterogeneous  medium  is  a  single  realization  of  an 
underlying  stochastic  process,  and  that  the  parameters  of  this  process 
may  be  approximated  using  the  statistics  of  the  conductivity  field. 
Also  from  these  statistics,  the  dispersive  capability  of  a  single  aqui¬ 
fer  may  be  determined  by  ensemble  averaging  over  the  conceivable  aquifer 
realizations  (Gillham  et  al.,  1984).  These  stochastic  models  have  been 
used  to  show  that  ensemble  solute  spreading  is  generally  not  Fickian 
(Gelhar  et  al.,  1979;  Matheron  and  de  Marsily,  1980;  Smith  and  Schwartz, 
1980;  Gelhar  and  Axness,  1983).  A  disadvantage  of  the  stochastic 
approach  is  that,  since  a  real  aquifer  is  conceived  as  a  single  realiza¬ 
tion,  the  solute  plume  must  migrate  a  sufficient  distance  so  that  the 
ensemble  averaging  used  in  the  stochastic  model  is  interpretable 
(Gillham  et  al.,  1984). 

Many  stochastic  models  that  have  been  developed  neglect  the  mecha¬ 
nisms  of  local  mechanical  dispersion  and  diffusion,  so  that  spreading  is 
strictly  a  consequence  of  heterogeneous  advection  (Warren  and  Sklba, 
1964;  Mercado,  1967;  Schwartz,  1977;  Smith  and  Schwartz,  1980).  A 
heterogeneous  advection  model,  when  applied  to  a  single  aquifer  realiza¬ 
tion,  predicts  local  Irregularities  in  the  concentration  distribution 
which  will  persist  at  a  macroscopic  scale  (Sudicky,  1983).  Sudicky 
(1983)  noted  that  in  a  field  experiment  where  detailed  measurements  of 
local  concentration  were  obtained,  this  phenomenon  was  not  observed,  and 
in  fact,  local  irregularities  observed  at  early  times  were  smoothed  out 


at  later  times.  Sudicky  (1983)  proposed  that  the  smoothness  of  the 
observed  macroscopic  concentration  patterns  was  the  result  of  transverse 
molecular  diffusion,  whereby  solute  which  advects  rapidly  in  the  high 
permeability  zones  diffuses  into  the  less  permeable  zones.  Sudicky 
(1983)  and  Gtlven  et  al.  (1984)  discussed  the  impact  of  diffusion  on 
solute  transport  in  a  layered  system,  and  compared  the  so-called  advec- 
tion/dlf fusion  model  with  stochastic  models.  It  was  shown  that  the 
deterministic  advection/dif fusion  model  provided  results  which  were 
equivalent  to  results  obtained  from  the  stochastic  theory  developed  by 
Gelhar  et  al.  (1979).  Sudicky's  (1983)  and  Gtlven  et  al.'s  (1984)  stu¬ 
dies  were  limited  to  the  analysis  of  conservative  solute  transport  in  a 
stratified  medium,  where  each  stratum  had  a  different  hydraulic  conduc¬ 
tivity,  and  therefore  a  different,  though  steady,  groundwater  flow 
velocity. 

The  importance  of  diffusion  into  spehrical  zones  of  low  permeabil¬ 
ity  was  demonstrated  in  column  experiments  conducted  by  Rao  et  al. 
(1980)  and  Nkedi-Kizza  et  al.  (1982).  Goltz  and  Roberts  (1984),  using 
numerical  simulations,  demonstrated  how  a  solute  transport  model  which 
neglected  the  mechanism  of  diffusion  could  significantly  underpredict 
long-time  contaminant  concentrations. 

The  effect  on  solute  transport  of  diffusion  into  low  permeability 
zones  has  been  discussed  in  the  chemical  engineering  literature  over  at 
least  the  past  thirty  years  (Deisler  and  Wilhelm,  1953;  Vermuelen,  1953; 
Rosen,  1954).  More  recently,  these  diffusion  concepts  have  been  applied 
to  study  contaminant  transport  by  groundwater  (van  Genuchten  and 
Wlerenga,  1976;  Rao  et  al.,  1980;  Nkedi-Kizza  et  al.,  1982;  Valocchi, 
1985a;  van  Genuchten,  1985;  Goltz  and  Roberts,  1986;  Crittenden  et  al., 
1986;  Miller  and  Weber,  1986).  These  studies,  however,  have  all  pre¬ 
sumed  one-dimensional  solute  transport.  Since  the  groundwater  environ¬ 
ment  is  three-dimensional,  there  is  value  in  deriving  three-dimensional 
formulations  of  transport  models  which  incorporate  the  diffusion  mecha¬ 
nism. 

SCOPE  OF  THIS  INVESTIGATION 

This  work  was  undertaken  in  response  to  the  apparent  importance  of 
diffusion  in  affecting  the  groundwater  transport  of  contaminants. 


Previous  models  which  have  incorporated  the  diffusion  mechanism  have 
assumed  one-dimensional  transport,  whereas  this  study  will  present  a 
three-dimensional  formulation.  The  specific  objectives  of  this  investi¬ 
gation  can  be  summarized  as  follows: 

1.  Develop  and  test  a  transport  model  which  incorporates  the 
mechanisms  of  advection,  dispersion,  linear  reversible 
sorption,  and  diffusion  in  a  three-dimensional,  infinite 
medium.  This  study  will  be  limited  to  the  special  case 
where  advection  is  due  to  a  single,  steady,  groundwater 
flow  velocity. 

2.  Compare  how  simulations  of  such  a  model  differ  from  simu¬ 
lations  of  the  advection/dispersion  model,  which  tradi¬ 
tionally  has  been  used  to  describe  contaminant  transport. 

To  facilitate  this  comparison,  develop  and  apply  methods 
for  obtaining  spatial  and  temporal  moments  of  the  concen¬ 
tration  distributions  simulated  using  the  different 
models. 

3.  Compare  and  contrast  simulations  of  the  different  diffu¬ 
sion  model  formulations. 

4.  Apply  the  model  to  the  data  set  obtained  in  an  extensive, 
high-resolution  field  experiment.  Assess  whether  char¬ 
acteristics  of  the  spatial  and  temporal  data  are  explain¬ 
able  using  models  which  incorporate  diffusion.  In 
addition,  examine  if  alternate  models  which  assume  other 
mechanisms  (e.g. ,  nonlinear,  hysteretlc  sorption)  can 
explain  experimental  observations.  This  will  require  an 
analysis  of  the  Bpatial  moment  behavior  of  these  alter¬ 
nate  models. 


CHAPTER  2 
MODEL  FORMULATION 


This  chapter  reviews  the  so-called  two-region  models,  which  couple 
an  expression  describing  advective/dispersive  solute  transport  in 
regions  of  mobile  water  with  an  expression  describing  mass  transfer  into 
regions  of  immobile  water. 

A  three-dimensional  formulation  for  such  two-region  models  will  be 
presented  and  an  analytical  solution  derived.  The  solution  technique 
Involves  consecutively  applying  the  Laplace  transform  in  time,  and  the 
Fourier  transform  in  space  to  the  coupled  set  of  partial  differential 
equations  that  mathematically  describe  the  two-region  models. 

REVIEW  OF  EXISTING  MODELS 

Transport  of  hydrophobic  organic  chemicals  in  groundwater  tradi¬ 
tionally  has  been  described  using  the  homogeneous  advective/dispersive 
transport  equation  with  a  sink  term  to  account  for  sorption  of  the 
organic  solute  onto  the  soil  matrix.  This  sorption  term  is  often  devel¬ 
oped  assuming  local  equilibrium  and  a  linear,  reversible  equilibrium 
relationship  between  the  quantity  of  the  organic  compound  in  the  sorbed 
and  solution  phases.  Several  investigators  have  found,  in  laboratory 
column  studies,  that  the  nearly  symmetric,  sigmoid  forms  of  breakthrough 
responses  predicted  using  models  making  these  simplifying  assumptions, 
do  not  agree  with  experimental  observations  (van  Genuchten  and  Wierenga, 
1976;  Rao  et  al.,  1979;  Reynolds  et  al.,  1982;  De  Smedt  and  Wierenga, 
1984).  Frequently,  experimental  breakthrough  responses  exhibit  highly 
asymmetric  or  nonsigmoid  profiles,  commonly  termed  tailing.  Tailing  may 
be  attributable  to  the  Blow  diffusion  of  solute  into  zones  of  immobile 
water.  It  has  been  hypothesized  that  these  zones  result  from  soil 
aggregation,  slow  flow,  or  unsaturated  flow  (van  Genuchten  and  Wierenga, 
1976;  Rao  et  al.,  1980;  Nkedi-Klzza  et  al.,  1982;  De  Smedt  and  Wierenga, 
1984). 

Various  models  have  been  proposed  to  describe  the  exchange  of 
solute  between  mobile  and  immobile  zones.  The  simplest  of  these,  the 
first-order  rate  model,  assumes  completely  mixed  zones  of  immobile 
water,  with  a  first-order  rate  expression  describing  diffusional 


transfer  of  solute  between  the  mobile  and  immobile  regions  (Coats  and 
Smith,  1964;  van  Genuchten  and  Wierenga,  1976).  These  models  couple  the 
advection/dispersion  equation  with  a  first-order  rate  expression. 
First-order  rate  models  have  successfully  simulated  the  observed  tailing 
In  laboratory  column  solute  transport  experiments  (van  Genuchten  and 
Wierenga,  1977;  van  Genuchten  et  al.,  1977;  De  Smedt  and  Wierenga, 
1979a;  Nkedi-Kizza  et  al.,  1984).  Analytical  solutions  for  this  type  of 
model  have  been  derived  for  different  initial  and  boundary  conditions 
applicable  to  finite  and  semi-infinite  columns  (Llndstrom  and  Narasimhan, 
1973;  van  Genuchten  and  Wierenga,  1976;  De  Smedt  and  Wierenga,  1979b). 

More  complex  models  have  been  developed  to  describe  the  transfer  of 
solute  within  immobile  regions  by  Flck's  second  law  of  diffusion.  These 
models,  which  couple  the  advection/dispersion  equation  with  an  expres¬ 
sion  to  describe  diffusion,  will  be  referred  to  as  diffusion  models. 
Diffusion  models  have  also  been  successfully  used  to  simulate  the  ob¬ 
served  tailing  in  laboratory  column  solute  transport  experiments  (Rao  et 
al.,  1980;  Nkedi-Kizza  et  al.,  1982).  These  models  assume  a  geometry 
for  the  immobile  region.  One-dimensional  analytical  solutions  to  diffu¬ 
sion  models  have  been  derived,  for  semi-infinite  boundary  conditions, 
assuming  spherical  (Pellett,  1966;  Rasmuson  and  Neretnleks,  1980), 
rectangular  (Sudicky  and  Frind,  1982),  and  cylindrical  (Pellett,  1966; 
van  Genuchten,  1985)  immobile  region  geometries.  Van  Genuchten  (1985) 
summarizes  solutions  for  different  immobile  region  geometries. 

Recent  research  has  begun  to  focus  on  solute  transport  in  "semi- 
controlled"  field  settings  (Leland  and  Hillel,  1982;  Roberts  et  al., 
1982;  Sudicky  et  al.,  1983;  Mackay  et  al.,  1986).  Solute  pulses  have 
been  Introduced  into  groundwater  under  conditions  corresponding  to  an 
infinite  medium,  in  which  the  medium  is  effectively  unbounded  upgradlent 
from  the  point  of  introduction  as  well  as  downgradient  from  the  polnt(s) 
of  observation.  To  permit  proper  analysis  of  the  data  from  the  perspec¬ 
tive  of  dispersion  and  diffusion  phenomena,  solutions  to  the  transport 
equation  under  the  pertinent  boundary  conditions  in  a  three-dimensional 
infinite  medium  are  required.  Although  one-,  two-,  and  three- 
dimensional  solutions  to  the  advection/dispersion  equation  in  an  infi¬ 
nite  porous  medium  are  available  (Carslaw  and  Jaeger,  1959;  Bear,  1972; 
Hunt,  1978),  of  the  two-region  models,  only  the  first-order  rate  model 
has  been  solved  analytically  for  infinite,  multi-dimensional  conditions 


(Carnahan  and  Remer,  1984).  Blbby  (1979)  combined  a  two-dimensional 
finite  element  model  with  an  analytical  expression  describing  solute 
diffusion  into  layers,  to  simulate  chloride  movement  in  a  chalk  aquifer. 
However,  to  date,  no  multi-dimensional  analytical  solutions  to  diffusion 
models  have  been  presented.  Such  multi-dimensional  solutions,  with 
infinite  boundary  conditions,  are  presented  in  this  chapter. 


ADVECTION/DISPERSION  MODEL  SOLUTION 

Sorbing  solute  transport  through  a  porous  medium  has  often  been 
described  using  the  advective/dispersive  transport  equation  with  a 
sorption  term  (Bear,  1972): 

3C(x , y ,  z ,  t )  x  _  ,  9  C  ,  _  ,  9 _ C  ,  _  i  9  C  9C  p  9S  ,  ~  . . 

9t  DX  3x2  y  ay2  az2  o  9x  ~  0  3t  (2"l> 

where  C(x,y,z,t)  represents  the  aqueous  solute  concentration,  S(x,y,z,t) 
is  the  sorbed  solute  concentration,  0  is  the  aquifer  porosity,  p  the 
bulk  soil  density,  D^,  Dy,  and  Dz  are  the  principal  components  of  the 
dispersion  tensor  in  the  x-,  y-,  and  z-directions,  respectively,  and  vQ 
is  the  average  pore  water  velocity.  Equation  2-1  implicitly  assumes 
steady,  uniform  flow  in  a  homogeneous,  isotropic  porous  medium.  If 
linear,  reversible,  equilibrium  sorption  is  assumed,  sorbed  and  aqueous 
solute  concentrations  may  be  related  using  the  concept  of  a  partition  or 
distribution  coefficient,  Kd,  such  that: 

S  -  KdC 

Uith  these  assumptions,  a  dimensionless  retardation  factor,  R,  can  be 
defined: 


»  ■  1  +  f  Kd 


so  Eq.  2-1  can  be  rewritten  as  (Bear,  1972): 
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Equation  2-2  has  been  solved  for  the  following  initial/boundary  condi¬ 
tions,  representing  an  instantaneous  point  source  in  an  infinite  medium: 
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C(x,y,z,0)  -  6(x)6(y)6(z) 


(2-3a) 


C(±»,y,z,t)  -  C(x,±»|2,t)  -  C(x,y,±«,t)  -  0 


(2-3b) 


The  solution  is  (Carslaw  and  Jaeger,  1959;  Crank,  1975;  Hunt,  1978): 
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FIRST-ORDER  RATE  MODEL  SOLUTION 


He  will  begin  our  discussion  of  the  two-region  models  by  deriving 
the  solution  to  the  first-order  rate  model  in  a  three-dimensional,  infi¬ 
nite  medium.  This  derivation  will  serve  as  a  guide  for  the  solution  of 
the  slightly  more  complex  diffusion  models. 

In  three  dimensions,  the  first-order  rate  model  for  sorbing  solute 
transport  in  a  porous  medium  with  immobile  water  zones  may  be  written 
(van  Genuchten  and  Wlerenga,  1976): 


3Cm(x,y,z,t) 

3t 


32C  D 

_ ®  +  -EL 

„  2  R 
3x  m 


jn  ^jn 

R  3x 
m 


ei»Rim  3Ci» 


6mRm 


3t 


(2-5) 


3C.  (x,y,z,t)  , 

iT2 - -  FT"  (Cm  ’  Cim> 

3t  0imRim  B  iB 


(2-6) 


These  equations  assume  groundwater  flow  in  the  positive  x  direction,  and 
that  D^,  Dmy,  and  Dmz  represent  the  mobile  zone  dispersion  coefficients 
in  the  x-,  y-,  and  z-directions.  Cm  and  represent  solute  concentra¬ 
tions  in  the  mobile  and  immobile  regions,  respectively,  and  6,,,  and 
are  mobile  and  immobile  region  water  content,  such  that  0  -  0m  +  0im* 
Equation  2-5  is  the  advection/disperslon  equation  with  a  sink  term  to 
describe  the  mass  transfer  of  solute  from  the  mobile  to  the  Immobile 
water  region.  Sorption  onto  the  solids  is  assumed  to  be  linear  and 
reversible,  with  the  effect  of  sorption  incorporated  into  Rm  and  Rim, 
the  retardation  factors  for  the  mobile  and  immobile  regions.  Following 
van  Genuchten  and  Wlerenga  (1976),  define: 
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and 


where  f  is  the  fraction  of  sorption  sites  adjacent  to  regions  of  mobile 
water.  Equation  2-6  is  the  first-order  rate  expression  describing 
solute  transfer  between  the  mobile  and  immobile  regions,  where  a'  is  the 
first-order  rate  constant. 

For  an  instantaneous  point  source  of  solute  in  an  infinite  porous 
medium,  the  following  initial/boundary  conditions  apply  to  2-5  and  2-6: 


jji 

Cn(x,y,z,0)  -  -g-g-  6(x)«(y)6(z) 


(2-7a) 


clm(x,y,z,°)  -  0  (2-7b) 

C  (±»,y,z,t)  -  C  (x,4«,z,t)  -  C  (x,y,±«,t)  -  0  (2-7c) 
n  in  id 

Equation  2-7a  represents  the  initial  solute  distribution  in  the  mobile 
zone,  (2-7b)  states  that  there  is  initially  no  solute  in  the  Immobile 
region,  and  (2-7c)  gives  the  boundary  conditions  for  an  infinite  medium. 

The  method  to  solve  (2-5)  and  (2-6)  simultaneously,  subject  to 
lnltlal/boundary  conditions  (2-7),  is  given  in  Appendix  A.  The  solution 
is 


i .  »  t 

Cm(x,y,*,t)  -  exp(T2^-)G(x,y,z,t)  +  -r-y  /  H(t,r)G(x,y,z,T>dT  (2-8) 
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®nd  IjJ  }  is  a  modified  Bessel  function.  This  solution  is  equivalent  to 
the  solution  presented  by  Carnahan  and  Remer  (1984). 

The  principle  of  superposition  can  now  be  used  to  integrate  the 
point  source  response  to  obtain  responses  for  other  initial  conditions. 
Carnahan  and  Remer  (1984)  provided  solutions  for  various  initial  condi¬ 
tion  geometries. 


DIFFUSION  MODEL  SOLUTION 

For  diffusion  models,  Eq.  2-5,  the  advection/dispersion  equation, 
is  still  applicable.  However,  since  these  models  allow  for  concentra¬ 
tion  gradients  within  the  Immobile  regions,  the  dependent  variable  C£m 
in  (2-5)  now  represents  a  volume-averaged  solute  concentration  within 
the  immobile  zone.  Of  course,  in  the  first-order  rate  model,  where  the 
Immobile  zone  is  assumed  to  be  perfectly  mixed,  the  solute  concentration 
throughout  the  immobile  zone  equals  the  average  concentration. 


Spherical  Geometry 

Equation  2-5  makes  no  assumption  regarding  the  immobile  zone  geome¬ 
try.  Assuming  spherical  Immobile  zones  of  radius  b,  the  volume  averaged 
concentration,  C^m,  can  be  defined  as  follows: 


3  k  2 

:i®  "  T  f  r  C  (r,x,y,z,t)dr 
b  0 


Fick's  law  of  diffusion  within  a  sphere  is  written: 
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with  the  boundary  conditions: 

Ca(0,x,y,z,t)  *  •  (2-lla) 

Ca(b,x,y,z,t)  -  (^(x.y.z.t)  (2-llb) 

The  method  to  solve  (2-5),  (2-9),  and  (2-10)  simultaneously,  for  an 
instantaneous  point  source,  subject  to  initial/boundary  conditions  (2-7) 
and  (2-11),  is  presented  in  Appendix  B.  The  solution  is: 
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Superposing  this  solution  will  provide  responses  for  other  initial  con¬ 
ditions.  It  should  be  noted  that  solving  for  Cm  in  Eq.  2-12  requires 
the  evaluation  of  an  infinite  integral.  The  integrand  of  this  Integral 
is  the  product  of  an  exponentially  decaying  function  and  a  sinusoidally 
oscillating  function.  The  Gaussian  quadrature  methods  described  by 
Rasmuson  and  Neretnieks  (1981)  and  van  Genuchten  et  al.  (1984)  may  be 
used  to  evaluate  the  integral  numerically. 

Cylindrical  Geometry 


For  cylindrical  Immobile  zones  of  radius  b,  the  volume  averaged 
concentration  is: 


(2-13) 


\m  "  72  /  rCa(r,x,y,*,t)dr 
b  0 


and  Fick's  law  of  diffusion  in  a  cylinder  may  be  written: 
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with  the  boundary  conditions  (2-lla)  and  (2-llb).  The  analysis  proceeds 
analogously  to  the  spherical  geometry  case,  with  the  three-dimensional 
point  source  response  the  same  as  (2-12),  but  in  the  case  of  cylindrical 


geometry: 
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Rectangular  Geometry 


In  the  case  of  a  rectangular  Immobile  zone,  it  is  physically  more  re¬ 
alistic  to  consider  advectlve/dispersive  transport  in  a  two-dimensional 
mobile  zone  with  diffusion  into  rectangular  immobile  zone  layers  of 
half-width  b.  The  set  of  equations  which  describes  such  a  system  is: 
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with  the  lnltial/boundary  conditions: 
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Using  the  methods  given  in  Appendix  B,  these  equations  can  be 
solved  to  give  the  following  result: 
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where: 
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Kq  is  a  modified  Bessel  function  of  the  second  kind,  and  Zp  and  Zm  are 
as  defined  in  Eq.  2-12,  though  for  rectangular  geometry: 
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MODEL  TESTING 


In  the  limiting  case,  as  the  volume  of  water  in  the  immobile  region 
becomes  small,  and  the  diffusion  rate  within  the  immobile  region  becomes 
large,  the  diffusion  model  solution  should  approach  the  solution  of  the 
advection/dispersion  model.  To  verify  this,  a  test  situation,  which  is 
depicted  in  Figure  2.1,  was  devised.  The  concentration  response  to  a 
rectangular  prism  initial  solute  distribution  of  half  length,  width,  and 
depth  L,  M,  N,  respectively,  may  be  calculated  at  a  sampling  well  using 
the  advection/dispersion  and  spherical  diffusion  models.  The  solution 
to  the  three-dimensional  advection/dispersion  model  for  a  rectangular 
prism  source  is  well  known  (Carslaw  and  Jaeger,  1959;  Hunt,  1978).  This 
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Figure  2.1.  Initial  condition  for  three-dimensional  model  testing. 
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solution,  for  realistic  parameter  values,  is  plotted  in  Figure  2.2,  as  a 
breakthrough  response  curve  at  the  sampling  well.  Also  shown  in  Figure 
2.2  is  the  solution  to  the  spherical  diffusion  model  with  very  low 
porosity  in  the  immobile  region  »  0.00038]  and  relatively  fast 
diffusion  [Dg/(R£mb^)  -  0.238  d-*]  within  the  immobile  region.  This 
solution  was  obtained,  for  a  rectangular  prism  initial  distribution,  by 
numerically  superposing  the  point  source  response  (Eq.  2-12)  over  the 
length,  width,  and  depth  of  the  initial  solute  distribution.  As  ex¬ 
pected,  the  solutions  are  identical,  because  the  fraction  of  porosity  in 
the  immobile  region  was  chosen  to  be  so  small  as  to  be  practically 
negligible. 

Another  means  of  testing  the  diffusion  model  solution  makes  use  of 
the  concept  of  approximate  equivalence  between  the  first-order  rate 
model  and  the  diffusion  model.  Van  Genuchten  (1985)  showed  that  the 
concentration  responses  for  the  two  models  would  be  approximately  the 
same  if  the  first-order  rate  parameter  and  the  spherical  diffusion  rate 
parameter  were  related  by  the  expression: 


o’  -  22.68 


D’0. 
e  1m 


(2-19) 


The  solution  to  the  first-order  rate  model  for  a  rectangular  prism 
source  may  be  obtained  by  superposing  the  point  source  response 
(Eq.  2-8)  over  the  length,  width,  and  depth  of  the  initial  solute  dis¬ 
tribution.  This  solution  is  Eq.  2-8  with: 
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Figure  2.3  shows  the  first-order  rate  model  solution  (Eqs.  2-8  and  2-20) 
for  some  realistic  parameter  values. 
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Figure  2.3.  Comparison  of  the  first-order  rate  and  spherical  diffusion 
model  solutions:  $  ■  0.90,v_  ■  0.091  m/d,  0  ■  0.38,  Dn  * 
0.02  mz/d,  D  -  0.0016  mZ/d,  D  -  0.6  x  10"  5  mz/d,  l  - 

5.0  m,  m  •  0.0  m,  n  »  -0.40  m,  L  »  1.5  m,  M  *  3.0  m,  N  * 
0.8  m,  -  2.78,  and  ■  5.00. 


Using  Eq.  2-19,  an  equivalent  spherical  diffusion  rate  parameter 
may  be  obtained  for  use  in  the  spherical  diffusion  model.  As  can  be 
seen  in  Figure  2.3,  the  solutions  of  the  first-order  rate  model  and  the 
spherical  diffusion  model  are  similar  in  form,  but  differ  slightly  in 
detailed  shape.  In  Chapter  4,  the  concept  of  model  equivalence  will  be 
discussed  in  more  detail. 
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CHAPTER  3 
MOMENT  ANALYSIS 


In  the  preceding  chapter,  models  were  presented  that  describe  solute 
transport  by  integrating  either  a  diffusion  expression  or  first-order 
rate  expression  into  the  three-dimensional  advective/dispersive  equa¬ 
tion.  A  convenient  means  of  quantitatively  studying  the  solute  plume 
behavior  predicted  using  such  two-region  or  physical  nonequilibrium 
models  is  to  examine  the  moments  in  space  and  time  of  the  models'  simu¬ 
lated  concentration  distributions.  In  this  chapter,  a  three-dimensional 
form  of  Arls'  method  of  moments  is  presented,  and  then  used  to  derive 
temporal  moments  associated  with  the  mobile  and  Immobile  regions.  In 
addition,  a  one-  and  three-dimensional  spatial  analog  to  Arls'  method  is 
developed,  and  used  to  examine  spatial  moment  behavior  in  both  the  mobile 
and  immobile  regions.  The  moment  analysis  is  extended  to  assess  the 
effect  of  model  dimensionality  on  the  form  of  the  moment  expressions. 


TEMPORAL  MOMENT  EQUATIONS 
Mobile  Region 

Based  on  the  definitions  of  the  jth  absolute  temporal  moment  of  a 
solute  concentration  distribution,  (^(x.t): 

AO 

m.  »  /  t^C  (x,t)dt  (3-1) 

J.t  ® 

and  the  Laplace  transform  of  the  function  CD(x,t): 


l[Ca(x,t)]  -  Cm(x,s) 


C  (x,t)dt 

m 


(3-2) 


where  s  is  the  Laplace  transform  variable,  Aris  (1958)  showed  that: 


(-1)J 


lim  [ 
8+0 


d%|B(x,s) 


ds 


J 


(3-3) 


Equation  3-3,  referred  to  as  Aris'  method  of  moments,  is  quite  useful, 
since  it  allows  the  calculation  of  temporal  moments  without  having  to 
Invert  the  Laplace  transform.  Arls'  method  has  been  widely  used, 
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particularly  in  chemical  engineering  research,  to  analyze  the  temporal 
moments  of  concentration  responses  simulated  by  models  which  combine 
one-dimensional  advective/dispersive  transport  with  a  diffusion 
expression  (Kucera,  1965;  Schneider  and  Smith,  1968;  Wakao  and  Kaguei, 
1982;  Valocchi,  1985a). 

Extension  of  Aris*  method  to  two  and  three  dimensions  is  straight¬ 
forward,  though  as  far  as  can  be  ascertained,  has  not  been  utilized 
previously  to  describe  solute  transport.  In  three  dimensions,  Eq.  3-3 
can  be  rewritten: 


.  d%  (x,y,z,s) 

(-l)J  urn  [—5 - t - ] 

S+0  dsJ 


(3-4) 


Table  3.1  presents  one-  and  three-dimensional  solutions,  in  the 
Laplace  domain,  for  the  mobile  region  concentration  distributions  of  the 
local  equilibrium,  first-order  rate,  and  diffusion  models.  These 
Laplace  domain  solutions  are  obtained  from  the  derivations  given  in 
Appendices  A  and  B. 

Equation  3-3  can  be  applied  to  the  one-dimensional  expressions  for 
mobile  solute  concentration  listed  in  Table  3.1a.  It  is  convenient  to 
present  the  resulting  moments  in  normalized  form,  where  the  jth  normal¬ 
ized  moment  is  defined  as: 


(3-5) 


Results  for  the  zeroth,  first,  and  second  normalized  absolute  temporal 
moments  are  presented  in  Table  3.2.  The  moments  for  the  local  equilib¬ 
rium  and  diffusion  models  have  been  reported  previously  by  Kucera 
(1965),  who  used  initial/boundary  conditions  identical  to  those  used 
herein. 

Similarly,  Eq.  3-4  may  be  applied  to  the  three-dimensional  expres¬ 
sions  in  Table  3.1b.  The  temporal  moments  for  the  three-dimensional 
models  are  presented  in  Table  3.3. 

An  examination  of  Tables  3.2  and  3.3  reveals  several  important  fea¬ 
tures  of  the  moments.  In  one  dimension,  the  zeroth  absolute  temporal 
moment,  m0>t,  is  constant  for  all  three  models,  and  is  equal  to  M^/v. 
The  v  term  in  the  denominator  is  due  to  the  model  initial  condition, 
which  is  expressed  as  a  Dirac  pulse  of  the  solute  in  space  (units  of  [M] ) , 
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|  TABLE  3.1 

|  MOBILE  REGION  SOLUTE  CONCENTRATION  IN  THE  LAPLACE  DOMAIN 

3  FOR  VARIOUS  MODELS 
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TABLE  3.2 

ABSOLUTE  TEMPORAL  MOMENTS  FOR  1-D  MOBILE  SOLUTE  CONCENTRATION  RESPONSES 


Local  Equilibrium 
Model* 

First-Order 

Rate  Model 

Diffusion 
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and  v  *  1  layered  diffusion 

v  *  2  cylindrical  diffusion 

v  *  3  spherical  diffusion 

*From  Kucera  (1965). 


whereas  the  zeroth  moment  of  a  temporal  distribution  (one  dimensional 
concentration  versus  time)  has  units  of  [ML~^T].  Thus,  to  transform  the 
zeroth  moment  of  the  initial  condition  in  space  [M]  to  the  zeroth  moment 
of  a  one-dimensional  temporal  distribution  [ML~*T],  it  is  necessary  to 
multiply  by  a  factor  with  units  [L“*T] .  That  factor  is  the  constant  1/v. 

As  expected,  the  zeroth  moment  for  the  one-dimensional  models  is 
independent  of  both  the  diffusion  rate  and  model  type,  since  all  the 
mass  which  was  initially  put  into  the  one-dimensional  space  must  even¬ 
tually  flow  past  the  sampling  point.  Perhaps  less  obvious  is  the  fact 
that  the  zeroth  moment  is  independent  of  the  diffusion  rate  and  model 
type  for  the  three-dimensional  models  as  well.  However,  considering 
a  pathway  connecting  the  initial  solute  distribution  with  any  partic¬ 
ular  sampling  location,  it  can  be  seen  that,  although  diffusion  would 
affect  the  speed  with  which  solute  "particles"  reach  the  sampling  point, 


TABLE  3.3 

ABSOLUTE  TEMPORAL  MOMENTS  FOR  3-D  MOBILE  SOLUTE  CONCENTRATION  RESPONSES 
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the  total  amount  of  solute  sampled  would  be  dependent  only  on  the  total 
injected  mass  (M3),  the  sampler  location  (£,m,n),  and  the  hydrodynamic 
parameters  (v,  Dx,  Dy,  Dz). 

It  is  interesting  to  compare  the  temporal  first  moments  obtained 
from  the  one-dimensional  and  three-dimensional  models.  Considering  the 
first  moment  of  the  three-dimensional  models  along  the  line  of  advective 
transport  (m  ■  n  ■  0),  it  is  found: 


(3-6) 


for  the  local  equilibrium  model,  and 


(3-7) 


1 


“i,t  -  7  <1+6) 

for  the  first-order  rate  and  diffusion  models.  These  values  are  less 
than  those  of  the  one-dimensional  models  by  a  constant  2Dx/v^  for  the 
local  equilibrium  model,  and  2Dx(l+0)/v^  for  the  physical  nonequilibrium 
models.  This  effect,  of  essentially  delaying  the  arrival  of  the  center 
of  mass  of  the  one-dimensional  models,  is  due  to  the  fact  that  with 
the  three-dimensional  models,  solute  which  disperses  in  the  negative 
x-direction  may  eventually  disperse  in  the  y-  and  z-directions  as  well, 
thereby  never  passing  the  sampling  point  on  the  y  »  z  •  0  axis.  How¬ 
ever,  with  the  one-dimensional  models,  the  solute  which  has  dispersed  in 
the  negative  x-direction  will  nevertheless  eventually  pass  the  sampling 
point,  thereby  increasing  the  first  temporal  moment.  For  the  same 
reason,  the  second  temporal  moments  of  the  one-dimensional  models  are 
greater  than  those  of  the  three-dimensional  models  along  the  line 
y  *  z  ■  0. 

Immobile  Region 

The  methods  of  the  preceding  section  may  be  applied  to  determine 
the  temporal  moments  of  solute  associated  with  the  Immobile  regions. 

Before  commencing  the  analysis,  however,  some  explanation  is  required  to 
define  what  is  meant  by  the  concentration  distribution  associated  with 
the  immobile  region. 

With  regard  to  the  temporal  moments,  the  immobile  region  concentra¬ 
tion  response  would  be  obtained  by  sampling  the  Immobile  region  at  a 
point  in  space.  Conceptually,  Imagine  a  sampling  point  which  yields 
volume-averaged  solute  concentrations  from  a  region  of  immobile  water. 

To  insure  mass  balance,  it  is  necessary  to  multiply  moments  obtained 
from  the  immobile  region  concentration  distribution  by  a  weighting 
factor.  This  weighting  factor  is  required  because  solute  is  unevenly 
distributed  between  the  mobile  and  immobile  regions.  To  determine  the 
value  of  this  weighting  factor,  compare  the  total  solute  masses  associ¬ 
ated  with  the  mobile  and  immobile  regions  of  an  incremental  volume  of 
aquifer  (dV) .  The  total  mass  associated  with  the  mobile  region  (aqueous 
plus  sorbed)  is  and  the  total  mass  associated  with  the  immo¬ 

bile  region  (aqueous  plus  sorbed)  is  CimeimRimdV .  The  ratio  of  immobile 
to  mobile  solute  masses  is  therefore: 

t 

( 

22  ! 


Thus,  to  insure  proper  mass  balance,  moments  obtained  from  the  immobile 
region  concentration  distribution  must  be  multiplied  by  the  weighting 
factor  6*  Therefore,  the  jth  temporal  moment  of  the  immobile  region 
concentration  distribution  may  be  defined  as: 


nj,t  "  /  tj[eCiB(x,t))dt 


(3-8) 


Of  course,  for  the  local  equilibrium  model,  6*0,  and  these  moments  are 
identically  zero. 

The  immobile  region  analogs  to  Eqs.  3-3  and  3-4  are: 


d-X  (x,s) 

n j  t  -  8(-l)J  lim  [ - - ] 

J,t  8-0  d8J 

njft  -  e(-l)J  li®  [ - - ] 

s-0  dsJ 


(3-9a) 


(3-9b) 


Equations  3-9a  and  3-9b  require  solutions  for  the  immobile  concentra¬ 
tions  in  the  Laplace  domain,  for  the  various  models.  Table  3.4  lists 
these  solutions  in  one  and  three  dimensions.  The  zeroth  and  first 
moments,  obtained  by  applying  Eqs.  3-9a  and  3-9b  to  the  immobile  solute 
concentration  expressions  listed  in  Table  3.4,  are  presented  in  Table 
3.5.  As  with  the  mobile  concentration  response  moments,  it  is  conve¬ 
nient  to  express  moments  greater  than  the  zeroth  in  terms  of  normalized 
moments,  defined  as: 


J.t  n 


hil 


(3-10) 


Compare  Tables  3.2  and  3.5,  to  find  that: 


“o  t  +  no  t  "  -  Ml 

O ,  t  0,t  V  1 
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for  the  one-dimensional  physical  nonequilibrium  models.  This  is  again 
due  to  the  conversion  of  an  initial  condition  in  space  to  a  zeroth 
moment  in  time.  It  should  be  noted,  however,  that  the  conversion  factor 


TABLE  3.4 


IMMOBILE  REGION  SOLUTE  CONCENTRATION  IN  THE  LAPLACE  DOMAIN 
FOR  VARIOUS  MODELS 
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where  C  (x,s)  for  each  model  is  defined  in  Table  3.1. 
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ID 


for  the  total  solute,  mobile  plus  Immobile,  is  (l+B)/v,  whereas  for  the 
mobile  solute  alone,  the  conversion  factor  was  1/v.  This  is  because  the 
total  solute  mass  movement  is  slower  than  the  mobile  region  mass  move¬ 
ment  due  to  the  Influence  of  the  immobile  region. 

This  effect  can  be  seen  more  clearly  by  comparing  the  first  tem¬ 
poral  moments  of  the  mobile  and  immobile  solute  concentration  responses. 
The  first  moment  of  the  immobile  region  response  lags  the  first  moment 
of  the  mobile  by  a  constant:  1/a  for  the  first-order  rate  model,  and 
1/a  for  the  diffusion  models. 


TABLE  3.5 


ABSOLUTE  TEMPORAL  MOMENTS  FOR  IMMOBILE  SOLUTE  CONCENTRATION  RESPONSES 
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SPATIAL  MOMENT  EQUATIONS 
Mobile  Region 

Analogous  to  definitions  (3-1)  and  (3-2),  the  one-dimensional  jth 
absolute  spatial  moment  of  the  concentration  distribution,  (^(x.t),  is: 


m1  -  /  xJC  (x,t)dx 


(3-12) 


and  the  Fourier  transform  of  the  function,  CD(x,t),  is  defined  by: 


where  p  is  the  Fourier  transform  variable.  Using  the  following  property 
of  the  Fourier  transform  (Spiegel,  1968): 


*  .  .  .  d^C  (p,t) 

*[*^Cm(x,t)]  -  /  e”lpx  x3Co(x,t)dx  -  iJ  - - - 

—oo  dp  i 


(3-14) 


take  the  limit  as  p  0  to  find: 


d^C  (p,t) 


/  xjC  (x.t)dx  -  ijlim  - 2L 

-•  p+0  dp 


(3-15) 


and  then  using  Eq.  3-12  write: 


.  d3C (p,t) 

m  -  ij  lim  [ - - - ] 

J  p+0  dpJ 


(3-16) 


This  is  the  one-dimensional  spatial  analog  to  Aris'  method  of  tem¬ 
poral  moments.  Although  Eq.  3-16  is  a  well  known  property  of  Fourier 
transforms  (Bremermann,  1963;  Bracewell,  1978),  it  apparently  has  never 
been  used  previously  to  obtain  moments  of  spatial  concentration  distri¬ 
butions  for  solute  transport  models. 

The  extension  of  Eq.  3-16  to  two  and  three  dimensions  is  straight¬ 
forward.  In  three  dimensions,  the  absolute  spatial  moment  of  the  con¬ 
centration  distribution  (^(x.y.Zjt)  is: 


eo  oo  OB 

III  xJy  z*Cm(x,y,z,t)dxdydz 


(3-17) 


The  Fourier  transform  in  three  dimensions  is  defined  (Bracewell,  1978) 
by 

HCm(*»y.*tt>]  -  Cjj(p,q,u,t)  - 

IJ!  e"i(px+qy+uz)  cjXty#z,t)dxdyd2  (3-18) 


where  p,  q,  and  u  are  the  Fourier  transform  variables  in  the  x-,  y- ,  and 
z-directlons,  respectively.  With  these  definitions,  follow  the  one- 
dimensional  analysis  directly  to  find: 

-  ■  i»*h  lt.  iL  S  [^(!.-.p.t.)1| 

p+0  dpJ  dqk  du£ 


(3-19) 


As  with  Arts'  method  for  temporal  moments,  these  spatial  moment  equa¬ 
tions,  (3-16)  and  (3-19),  are  quite  useful,  since  they  allow  the  calcu¬ 
lation  of  moments  in  the  Fourier  domain,  thereby  eliminating  the  need 
for  complicated  inversions  of  the  transforms. 

Appendix  C  presents  the  details  of  deriving  the  absolute  spatial 
moments  for  the  local  equilibrium,  first-order  rate,  and  diffusion 
models.  The  method  followed  for  all  three  models  is  essentially  the 
same.  The  model  equations,  with  appropriate  Initial/boundary  conditions 
for  a  solute  pulse  in  an  infinite  medium,  are  Fourier  transformed.  The 

A 

equations  are  then  solved  in  the  Fourier  domain  for  Cm(p,q,u,t). 
Equation  3-19  is  then  applied  to  obtain  the  zeroth,  first,  and  second 
absolute  spatial  moments.  The  derivations  were  done  using  the  three- 
dimensional  formulae.  Conversion  of  the  three-dimensional  results  to 
one  dimension  is  trivial,  using  the  relationships: 

°0  "  “000  (3-20a) 
“l  “  “100  (3-20b) 
“2  *  “200  (3— 20c) 

Thus,  in  space,  model  dimensionality  has  little  effect  on  the  moments. 
However,  in  time,  model  dimensionality  plays  an  Important  role. 

Table  3.6  presents  the  spatial  moments  for  the  three  models,  where 


(3-21) 


defines  t-he  normalized  absolute  moment.  The  local  equilibrium  model 
results  are  well  known  (e.g.,  Freyberg,  1986).  Valocchi  (1985b),  using 
a  method  different  from  that  presented  here,  derived  the  moments  in 
space  for  the  first-order  rate  model.  The  diffusion  model  moments 
presented  in  Table  3.6  are  new. 


Immobile  Region 


The  methods  of  the  preceding  section  may  be  applied  to  determine 
the  spatial  moments  of  solute  concentration  associated  with  the  immobile 
region.  Analogous  to  Eq.  3-8,  which  defines  the  temporal  moments  of 
immobile  region  solute,  we  may  define  the  spatial  moments  for  the  immo¬ 
bile  region  as: 


TABLE  3.6 

SPATIAL  MOMENTS  FOR  THREE-DIMENSIONAL  MOBILE  SOLUTE  DISTRIBUTIONS 


% wp.wia-n w vm j* tn *w nm  ««ML*aiMja ?f?ji •*3ai}3*s*x,79Jt*si??~tL~P'*  -"y-^srrrm r?i_  !wpmt«wi 


! 


and 


°jk£  "  /  /  /  xjykz*[BCim(x,y,z,t)]dxdydz 

J  —00  —00  —oo 


(3-22) 


Again,  for  the  local  equilibrium  model,  6  ■  0,  and  the  immobile  region 
moments  are  identically  zero.  The  Immobile  region  analog  to  Eq.  3-19 
is: 


jk£ 


8  i 


j+k+£ 


dJ  fdk  rd  elm<P»9.u.*> 

lim  — r  1—  [ - t - ] } 

p+0  dp  dq  du* 

q+0 

u+0 


(3-23) 


Equation  3-23  may  be  applied  to  the  Fourier  transformed  immobile  region 
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concentration,  C£m(p,q,u,t) ,  to  obtain  the  immobile  region  moments. 
Details  of  deriving  the  immobile  region  spatial  moments  for  the  first- 
order  rate  and  diffusion  models  are  presented  in  Appendix  D.  Expres¬ 
sions  for  the  zeroth  and  first  immobile  region  moments  are  listed  in 
Table  3.7.  Valocchi  (1985b)  derived  the  first-order  rate  expressions 
shown  in  Table  3.7  by  a  different  method. 

Testing 

Evaluation  of  the  temporal  moments  for  all  three  models,  as  well  as 
evaluation  of  the  spatial  moments  for  the  local  equilibrium  and  first- 
order  rate  models,  is  straightforward;  and  in  fact  can  easily  be  accom¬ 
plished  using  a  hand-held  calculator.  However,  calculation  of  the 
diffusion  model  spatial  moments  requires  numerical  evaluation  of  an 
infinite  integral  using  the  Gaussian  quadrature  technique  of  Rasmuson 
and  Neretnieks  (1981)  and  van  Genuchten  et  al.  (1984);  as  well  as  evalu¬ 
ating  a  limit  as  the  parameter  e  approaches  zero.  Therefore,  it  is 
essential  to  assess  the  accuracy  of  the  numerical  evaluation,  as  well  as 
to  determine  the  appropriate  values  of  e  to  use. 

One  way  to  test  the  solution's  accuracy  is  to  compare  the  solution 
with  the  known  solution  of  a  limiting  case.  As  8  approaches  zero,  the 
diffusion  model' 8  spatial  moments  should  approach  the  known  spatial  mo¬ 
ments  of  the  local  equilibrium  model.  Figure  3.1  shows  the  relative  error 
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TABLE  3.7 


i^TS' 


and  Q2»  and  a  are  88  defined  In  Appendix  C.  A  and  B,  are  as  defined  in  Table  3.6. 
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Figure  3.1.  Relative  error  of  the  spherical  diffusion  model  mobile  re¬ 
gion  spatial  moment  calculation  (comparison  with  limiting 
case  of  equilibrium  transport):  v  *  0.1  m/d,  Dx  =  0.04 
m^/d,  De/b^  *  0.006  d“^,  8  •  0.0001,  and  t  *  100  d. 


versus  log(c)  for  the  spherical  diffusion  model's  zeroth,  first,  and 
second  spatial  moments  in  the  mobile  region,  where  relative  error  is 
defined  as: 


Relative  Error  - 


LE  moment  -  Diffusion  moment 


(3-24) 


LE  moment 

Since  the  parameter,  e,  has  units  of  [T-*]  (see  Table  3.6),  it  is 
Important  that  e  be  small  with  respect  to  the  other  rate  constants  in 
the  solution:  a,  a/8,  and  1/t.  Thus,  the  constraint  is  imposed  that: 


e  «  Min(a,  a/6,  1/t) 


(3-25) 


For  the  parameters  used  in  Figure  3.1, 


Min(a,  a/B,  1/t)  ■  0.01 


(3-26) 


From  Figure  3.1,  it  is  seen  that  the  relative  error  of  the  zeroth  moment 
solution  approaches  zero  for  log  e  <  (-4).  For  the  first  moment,  if 
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(-8)  <  log  e  <  (-4)  the  relative  error  is  small,  and  for  the  second 
moment  solution,  If  (-5)  <  log  e  <  (-4),  the  relative  error  is  small. 
For  the  zeroth,  first,  and  second  moment  solutions,  an  upper  bound  on  e 
may  be  defined  by  the  expression: 


e  <  0.01[Min(a,  a/6,  1/t)] 


(3-27) 


However,  to  evaluate  the  first  and  second  moments,  there  also  is  a 
lower  bound  on  e,  as  the  definite  integrals  in  the  first  and  second 
moment  formulae  approach  minus  infinity  as  the  lower  limit  of  integra¬ 
tion  approaches  zero.  Thus,  at  very  small  values  of  c,  numerical  errors 
in  the  quadrature  technique  used  to  evaluate  the  definite  integrals 
become  substantial.  This  is  not  the  case  for  the  zeroth  moment  Integral; 
consequently,  there  is  no  lower  bound  on  e  in  the  zeroth  moment  evalua¬ 
tion  . 

The  lower  bound  of  e  in  the  first  and  second  moment  solutions  is 
dependent  not  only  on  the  model  parameters,  but  also  on  machine  preci¬ 
sion.  In  practice,  to  be  confident  of  a  first  or  second  moment  value 
obtained  using  the  formulae  in  Table  3.6,  it  is  necessary  to  calculate 
the  solution  over  a  range  of  values  for  e,  thus  obtaining  a  plot  similar 
to  that  shown  in  Figure  3.1.  A  flat  region  will  usually  be  found,  where 
the  moment  will  not  vary  over  a  range  of  one  or  two  orders  of  magnitude 
for  e«  This  value  for  the  moment  is  correct  within  a  few  percent,  as 
can  be  seen  from  Figure  3.1. 

Another  test  is  to  compare  the  moments  calculated  using  the  formu¬ 
lae  in  Table  3.6  with  the  moments  calculated  by  actually  computing  the 
spatial  distribution  of  the  solute,  using  formulae  derived  in  Chapter  2, 
at  repeated  points  in  space,  and  then  using  a  numerical  quadrature  tech¬ 
nique  to  obtain  the  moments  of  the  distribution.  Figure  3.2  gives  a 
plot  of  the  relative  error  versus  log(e)  for  the  spherical  diffusion 
model,  where  the  relative  error  is  now  defined  as: 


Relative  error  - 


Moment  from  distribution  -  Moment  from  formula 
Moment  from  distribution 


(3-28) 


In  comparing  the  errors  of  the  moment  estimator  (Figures  3.1  and  3.2), 
it  is  apparent  that,  with  the  zeroth  and  first  moments,  an  acceptable 
accuracy  is  assured  over  many  orders  of  magnitude  of  e>  However,  the 
window  of  acceptable  accuracy  is  much  narrower  in  the  case  of  the  second 
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Figure  3.2.  Relative  error  of  the  spherical  diffusion  model  mobile 
region  spatial  moment  calculation  (comparison  with 
numerically  evaluated  moments):  v  *  0.1  m/d,  Dx  »  0.004 
m2/d,  D  /b2  -  0.006  d”1,  6  ■  0.25,  and  t  ■  100  d. 


moment,  being  confined  to  1  to  2  orders  of  magnitude  in  e  in  the  example 


shown. 


A  similar  plot  can  be  obtained  for  the  zeroth  and  first  immobile 
region  moments,  by  comparing  the  formulae  for  the  spherical  diffusion 
model  (Table  3.7)  with  results  obtained  by  numerically  integrating  the 
Immobile  region  concentration  distribution.  The  immobile  region  concen¬ 
tration  distribution  was  obtained  by  using  the  IMSL  subroutine  FLINV 
(IMSL,  1982)  to  numerically  invert  the  Laplace  domain  solution  for  the 
Immobile  region  concentration  (Table  3.4),  at  repeated  points  in  space. 
Figure  3.3  gives  the  relative  error  versus  log(e)  for  the  zeroth  and 
first  immobile  region  moments  of  the  spherical  diffusion  model. 

For  the  parameter  values  used  in  Figures  3.2  and  3.3: 


Min(a,  a/0,  1/t)  ■  0.01 


(3-29) 


Again, 


e  <  0.01 [Min(a, a/6,  1/t)] 


(3-30) 
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Figure  3.3.  Relative  error  of  the  spherical  diffusion  model  immobile 
region  spatial  moment  calculation  (comparison  with  numeri¬ 
cally  evaluated  moments):  v  *  0.1  m/d,  Dx  ■  0.004  m2/d, 
De/b2  «  0.006  d"*,  6  *  0.25,  and  t  *  100  d. 

defines  an  upper  bound  on  e.  Note,  too,  the  sensitivity  of  the  first 
and  second  moment  solutions  to  excessively  small  values  of  e* 

In  summary,  the  numerical  evaluation  of  the  zeroth  moment  formulae, 
shown  in  Tables  3.6  and  3.7  for  the  geometrical  diffusion  models,  is 
straightforward,  as  long  as  e  is  two  orders  of  magnitude  less  than 
Min[a,  a/6,  1/t].  Evaluation  of  the  first  and  second  moments  is  not  as 
straightforward,  since  e  has  a  lower  bound  criterion  which  must  also  be 
satisfied,  and  which  is  not  easily  defined.  When  evaluating  the  first, 
and  particularly  the  second  moment  formulae,  it  is  necessary  to  deter¬ 
mine  how  the  result  varies  over  a  range  of  e,  before  accepting  a  value 
for  the  moment.  In  some  cases,  there  is  a  possibility  that  no  value  of 
c  can  be  found  which  meets  both  lower  and  upper  bound  criteria.  This  is 
especially  true  when  evaluating  the  second  moment,  where,  judging  from 
Figures  3.1  and  3.2,  the  lower  bound  of  e  can  be  orders  of  magnitude 
greater  than  the  lower  bound  of  the  first  moment  formula. 
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TEMPORAL  VERSUS  SPATIAL  DATA 


Using  the  relationships  derived  above,  which  permit  determination 
of  temporal  and  spatial  moments  given  parameters  for  the  three  models 
under  discussion,  it  is  possible  to  find  how  certain  effective  parame¬ 
ters  behave  with  respect  to  model  type,  and  with  respect  to  type  of 
measurement  (in  space  or  time).  In  particular,  an  effective  velocity 
(vgff)  and  an  effective  dispersion  coefficient  (Deff)  will  be  derived  in 
terms  of  model  parameters.  These  effective  parameters  are  local  equi¬ 
librium  model  equivalents  which  approximately  duplicate  the  concentra¬ 
tion  responses  of  the  physical  nonequilibrium  models.  Use  of  these 
effective  parameters  will  aid  in  understanding  how  the  general  behavior 
of  the  spatial  and  temporal  concentration  distributions  differ. 


Behavior  of  Mobile  and  Immobile  Temporal  Responses 

Effective  velocity  and  dispersion  coefficients  may  be  defined  in 
terms  of  the  temporal  moments  as  follows  (Valocchi,  1985a): 


and 


(3-31a) 


where 


M2,t  2 
2yl ,t  "eff 


v2,t  *  u2,t  ~  <vi,t^ 


(3-31b) 


is  the  central  second  moment . 

Thus,  using  the  expressions  for  the  mobile  region  temporal  moments 
in  Tables  3.2  and  3.3,  it  is  straightforward  to  determine  ve£f  and  Dejf 
in  terms  of  the  parameters  of  the  one-  and  three-dimensional  local  equi¬ 
librium,  first-order  rate,  and  diffusion  models. 

Table  3.8a  lists  the  effective  parameters  in  one  dimension  for  the 
three  models,  making  the  assumption  of  large  Peclet  number  (Pe  *  vfc/Dx 
»  1).  Table  3.8b  lists  the  effective  parameters  in  three  dimensions, 
along  the  line  y  *  z  *  0.  Due  to  the  assumption  of  large  Peclet  number 
in  the  one-dimensional  case,  the  effect  of  dispersion  in  the  negative 
x-direction  (which  was  discussed  earlier)  becomes  negligible,  so  that 
the  one-  and  three-dimensional  solutions  are  identical.  Of  course,  since 
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TABLE  3.8 

EFFECTIVE  VELOCITY  AND  DISPERSION  COEFFICIENTS  FROM  TEMPORAL  MOMENTS 


Local  Equilibrium  First-Order  Rate 

Model  Model 


Diffusion 

Model 


a.  One  Dimension  (at  high  values  of  Pe) 


°«  .  v28 

1+6  a(l+B)3 


v(v+2)D 


where  a  “ 


and  v  is  as  defined  in  Table  3.3. 


b.  Three  Dimensions  (along  the  line  y  ■  z  ■  0) 


Dx  .  v2B 
1+6  a(l+B)3 


,  v2B 
1+6  a(l+B)3 


,  v2B 
1+6  a(l+B)3 


the  effective  parameters  represent  local  equilibrium  model  equivalents, 
ve£j  and  for  the  local  equilibrium  model  are  v  and  Dx,  respectively. 
The  high  Peclet  number  expressions  for  vgf£  and  Deff  for  the  first-order 
rate  and  diffusion  models,  which  are  presented  in  Table  3.8a,  are  well 
known  in  one  dimension  (Baker,  1977;  De  Smedt  and  Wierenga,  1984; 
Valocchi,  1985a).  The  application  to  three  dimensions  is  new.  By 
applying  Eq.  3-31  to  the  moment  formulae  of  Table  3.3,  values  for  effec¬ 
tive  velocity  and  dispersion  coefficients  can  be  obtained  at  any  point 
in  space,  for  a  given  model. 

The  effective  velocity  of  the  immobile  region  response  (vneff)  may 
be  obtained  by  applying 


vnef  f 


(3-32) 


to  the  immobile  region  moment  formulae  shown  in  Table  3.5.  Again  con¬ 
sidering  high  Peclet  numbers  in  the  one-dimensional  case,  and  along  the 
line  y  ■  z  -  0  in  the  three-dimensional  case,  yields: 


'neff  *  £(1+6)  +  \ 


(3-33a) 


for  the  first-order  rate  model,  and 


neff  1(1+6)  +  l 
v  a 


(3-33b) 


for  the  diffusion  models.  Of  course,  the  concept  of  an  Immobile  region 
solute  distribution  velocity  is  not  applicable  to  the  local  equilibrium 
model.  As  expected,  the  effective  velocity  in  the  immobile  region  is 
less  than  that  in  the  mobile  region. 

In  a  forthcoming  section,  the  equations  derived  in  this  section 
will  be  used  to  compare  temporal  and  spatial  behavior. 


Behavior  of  Mobile  and  Immobile  Spatial  Distributions 


The  following  relationships  between  the  spatial  moments  and  the 
total  mass  associated  with  the  mobile  region  (MT),  the  center  of  mobile 
mass  location  (xc,  yc,  zc)  and  the  elements  of  the  spatial  covariance 
tensor  may  be  written  (Freyberg,  1986): 
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(3-34) 
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As  with  the  temporal  moment  discussion,  effective  velocity  and  disper¬ 
sion  coefficients  may  also  be  defined  as: 


i  I  »-k. »  ^  ♦  Ji  iJS  '-•mi  ««t_ 
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(3-35) 

(3-36) 


Thus,  starting  with  the  spatial  moment  relationships  of  Table  3.6,  it  is 
possible  to  derive  expressions  for  the  effective  velocity  and  dispersion 
coefficients  for  the  three  models  under  consideration.  In  practice, 
evaluating  the  effective  velocity  and  dispersion  coefficients  for  the 
diffusion  models  requires  the  numerical  evaluation  of  the  derivatives  in 
(3-35)  and  (3-36). 

Figure  3.4  shows  how  Hj.,  ve££,  and  Deffx  behave  as  a  function  of 
time  for  each  of  the  three  model  types.  Model  parameters  for  the  exam¬ 
ple  were  arbitrarily  chosen,  except  that  the  first-order  (a)  and  the 
diffusion  model  (a)  rate  parameters  were  selected  so  that 


a  •  a 


(3-37) 


The  use  of  this  equality  to  obtain  the  first-order  rate  parameter  will 
be  discussed  further  in  the  next  chapter. 

Figure  3.4a  compares  the  total  solute  mass  associated  with  the  mo¬ 
bile  region  (M^)  versus  time  for  each  of  the  models.  In  the  local  equi¬ 
librium  model,  the  total  mass  is  constant  over  time.  The  physical  non- 
equilibrium  models  show  a  loss  of  mass  with  time.  Initially,  all  the 
mass  is  associated  with  the  mobile  region.  With  time,  more  and  more 
solute  diffuses  into  the  immobile  region,  so  that  eventually,  the  mass 
in  the  mobile  region  equals  1/(1+B)  times  the  original  mass.  The  dif¬ 
ferences  between  the  first-order  rate  and  diffusion  models  do  not  appear 
to  be  significant. 

The  behavior  of  ve££  shown  in  Figure  3.4b  is  similar  to  the  total 
mass  behavior.  The  local  equilibrium  model  predicts  constant  velocity 
(ve££  -  v) ,  whereas  the  physical  nonequilibrium  models  predict  a  decel¬ 
eration  with  time,  going  from  veff  -  v  at  t  »  0  to  veff  -  v/(l+g)  at 
large  times. 

It  is  important  to  realize  that  the  simulated  decline  in  mass  and 
deceleration  of  the  mobile  region  plume  are  due  to  the  initial  condi¬ 
tion,  which  assumes  no  mass  associated  with  the  Immobile  region  at  t  *  0. 
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Figure  3.4.  Behavior  with  time  of  the  mobile  region  plume  concentration 
distribution  for  equilibrium,  first-order  rate,  and  diffu¬ 
sion  models:  v  -  0.5  m/d,  Dx  -  0.004  m2/d,  <*  -  1.5  d"1, 
and  g  ■  0.5.  a)  Relative  mass,  b)  Effective  velocity,  and 
c)  Effective  dispersion  coefficient. 
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If  this  assumption  is  inappropriate,  the  mobile  region  plume  mass  and 
velocity  behavior  would  differ  from  the  Figure  3.4  simulations.  The 
validity  of  this  initial  condition  assumption,  for  a  field  experiment, 
will  be  discussed  in  Chapter  5. 

The  behavior  of  Dej£x  with  time  is  shown  in  Figure  3.4c.  The 
effective  dispersion  coefficient  in  the  x-direction  is  constant  for  the 
local  equilibrium  model  (De£fx  •  Dx) ,  while  the  values  of  Deffx  calcu¬ 
lated  using  the  physical  nonequilibrium  models  increase  from  Deffx  *  Dx 
at  t  ■  0  to 


(3-38) 


(3-39) 


for  the  diffusion  models,  at  large  times. 

Expressions  similar  to  Eqs»  3-34  and  3-35  may  be  written  relating 
the  total  mass  associated  with  the  immobile  region  (NT),  center  of 
immobile  mass  location  (xnc,  ync,  *nc),  and  the  effective  velocity  for 
the  immobile  distribution  (vne£f),  to  the  immobile  region  moments: 

NT  "  0mRmIlOOO 


Figure  3.5a  shows  the  total  mobile  and  immobile  masses  over  time 
for  the  first-order  rate  model.  Since  total  mass  must  be  conserved,  the 
sum  of  the  mobile  and  immobile  masses  remains  constant. 

Figure  3.5b  plots  the  effective  velocity  of  the  solute  plumes  in  the 
mobile  and  immobile  regions.  The  figure  illustrates  the  deceleration  of 
the  mobile  solute  distribution  from  vefj  ■  v  at  early  times  to  veff  * 
v/(l+8)  at  later  times.  The  immobile  solute  distribution  has  an  effec¬ 
tive  velocity  (vne£f)  of  v/2  initially.  This  Initial  value  is  a  con¬ 
sequence  of  the  first-order  rate  expression,  which  controls  the  amount  of 
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solute  entering  and  leaving  the  immobile  region.  At  early  times,  there 
is  a  large  concentration  gradient  driving  the  solute  into  the  immobile 
region,  but  no  gradient  in  the  reverse  direction,  since  the  initial 
immobile  region  concentration  is  assumed  equal  to  0.  Since  only  times 
very  close  to  t  *  0  are  being  considered,  the  concentration  of  the 
mobile  plume  may  be  assumed  approximately  constant  (i.e.,  very  little 
solute  is  transferred  into  the  immobile  region).  Thus,  based  on  the 
first-order  rate  expression  which  controls  the  transfer  of  solute,  the 
(small)  amount  of  solute  which  is  transferred  into  the  immobile  region 
between  time  t  *  0  and  t  ■  At  (as  the  mobile  plume  moves  from  x  ■  0  to 
x  *  ax)  is  approximately  the  same  as  the  (small)  amount  of  solute  which 
is  transferred  between  t  *  At  and  t  ■  2At  (as  the  mobile  plume  moves 
from  x  ■  ax  to  x  ■  2ax).  Therefore,  in  the  time  it  takes  the  mobile 
plume  to  move  a  distance  2Ax,  the  center  of  mass  of  the  immobile  plume 
moves  from  x  -  0  (at  t  -  0)  to  x  -  Ax  (at  t  -  2At),  leading  to  the 
result.  Indicated  in  Figure  3.5b,  that  as  t  0,  vneff  -*•  v/2.  In  a 
mathematical  analysis  of  two-layer  flow  incorporating  first-order  mass 
transfer  between  the  layers,  Christodoulou  (1986)  obtained  this  same 
result . 

Eventually,  at  long  times,  the  mobile  and  immobile  region  plumes 
must  move  with  the  same  velocity: 

vef f  "  vnef f  *  ITg  (3-41) 

Depending  on  the  values  of  v  and  B,  the  immobile  region  plume  will 
either  accelerate  or  decelerate  from  its  initial  value  of  v/2  to  attain 
this  final  velocity.  It  should  be  noted  that  if 


I  1+B 


(3-42) 


(i.e.,  B  ■  1),  the  immobile  region  plume  moves  at  a  constant  velocity 
for  all  time.  Figure  3.6  illustrates  this  behavior  for  the  special  case 
of  Eq.  3-42. 

Figure  3.7  is  the  spherical  diffusion  model  analog  of  Figure  3.5. 
The  total  mobile  and  immobile  solute  masses  of  the  diffusion  model 
(Figure  3.7a)  qualitatively  behave  as  did  the  masses  simulated  using 
the  first-order  rate  model  in  Figure  3.5a.  The  effective  velocities  of 
the  mobile  and  immobile  region  plumes  (Figure  3.7b)  also  behave  similarly, 
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Figure  3.6.  Comparison  of  the  mobile  and  immobile  region  plume  effec¬ 
tive  velocity  behavior  with  time  simulated  by  the  first- 
order  rate  model  with  6  ■  1.0. 


except  that  the  initial  Immobile  plume  velocity  is  2v/3  in  the  case  of 
the  diffusion  model.  This  difference  is  due  to  the  different  model 
structure  describing  solute  transfer  into  and  within  the  Immobile 
region.  Again,  at  large  times,  both  the  mobile  and  Immobile  region 
plumes  attain  an  effective  velocity  of  v/(l+8).  If 
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(3-43) 


a 


(i.e.,  8  ■  0.50),  then  the  immobile  plume  center  of  mass  moves  at  the 
same  velocity  at  short  and  long  times.  Figure  3.7b  illustrates  that  for 
8  ■  0.50  the  short-  and  long-time  velocities  are  equal,  while  at  inter¬ 
mediate  times  the  velocity-time  relationship  exhibits  a  shallow  minimum. 
This  minimum  is  due  to  the  mathematical  formulation  of  the  diffusion 
model. 


.'i 


43 


I 


EFFECTIVE  VELOCITY  (m/d)  RELATIVE  MASS 


mjiiuuir  i  mw’.i'n  j  g  yggggTTTgggreryTgg 


The  expressions  for  Hj,  veff,  Deff,  NT>  and  vneff  derived  in  the 
preceding  sections  provide  a  convenient  tool  for  comparing  the  temporal 
and  spatial  behavior  of  the  concentration  distributions.  The  three- 
dimensional  models  can  be  assessed  by  comparing  their  spatial  behavior 
over  the  entire  distribution  with  the  temporal  behavior  along  the  line 
y  *  z  *  0.  This  is  equivalent  to  comparing  the  result  of  sampling  the 
entire  concentration  distribution  at  points  in  time,  with  the  break¬ 
through  responses  obtained  from  sampling  wells  along  the  centerline  of 
the  solute  trajectory. 

Zeroth  moment 

Mobile  region.  Examining  the  zeroth  temporal  moment  of  the  mobile 
distributions,  it  is  seen  that 

M3 

%  *  —  —  (3-44) 

4u  /D  D  i. 

y  z 

whereas  the  zeroth  spatial  moment  is  model  dependent.  The  zeroth  spa¬ 
tial  moment  for  the  local  equilibrium  model  remains  constant  at 

M,.  -  M$  (3-45) 

while  the  zeroth  spatial  moments  of  the  physical  nonequilibrium  models 
decline  from 

Mj  -  M$  (3-46) 


at  early  times  to 

M3 
“  l+e 

at  later  times. 


(3-47) 


Immobile  region.  For  the  physical  nonequilibrium  models,  the  tem¬ 
poral  zeroth  moments  are 


4*  /D  D  t 

y  * 


(3-48) 


while  the  spatial  moments  increase  from 


NT  -  0 


(3-49) 


at  early  times  to 


(3-50) 


at  later  times.  Another  significant  difference  is  that  the  sum  of  the 
mobile  and  immobile  zeroth  spatial  moments  remains  constant  over  time, 
whereas  the  sum  of  the  temporal  moments  is  inversely  proportional  to  the 
sampling  distance,  £. 


First  moment 


Mobile  region.  Before  comparing  effective  velocities,  which  are 
derived  from  the  spatial  or  temporal  first  moments  of  a  distribution,  it 
is  useful  to  recall  the  definition  of  the  effective  velocity  (vgff). 
The  effective  velocity  is  defined  as  the  velocity  parameter  which  would 
be  used  in  an  equilibrium  model  to  obtain  the  same  spatial  or  temporal 
first  moment  that  would  be  calculated  using  an  equivalent  nonequilibrium 
model.  Effective  velocity  does  not  represent  a  physical  quantity,  and 
as  this  section  will  emphasize,  the  effective  velocity  obtained  from  an 
analysis  of  the  spatial  distribution  is  both  quantitatively  and  qualita¬ 
tively  different  from  the  effective  velocity  calculated  from  a  temporal 
moment  analysis. 

Comparing  effective  velocities,  it  is  seen  that  for  the  physical 
nonequilibrium  models,  the  effective  velocity  obtained  using  the  tem¬ 
poral  moments  remains  constant  at 


vef f  "  ITg  (3~51) 

at  all  sampling  distances,  while  the  effective  velocity  obtained  using 
spatial  moments  declines  from 

vgff  -  v  (3-52) 


at  early  times  to 

vef f  “  TTf 


(3-53) 


at  later  times.  This  difference  can  be  understood  qualitatively  in  the 
following  way:  since  the  temporal  first  moment  reflects  behavior  aver¬ 
aged  over  all  time,  it  provides  an  equilibrium  value  for  the  effective 
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velocity,  whereas  the  spatial  first  moment  estimate  provides  a  "snapshot" 
of  the  instantaneous  velocity  with  which  the  solute  distribution  is 
moving,  at  particular  points  in  time. 

This  difference  in  spatial  and  temporal  behavior  has  interesting 
repercussions  in  light  of  the  concept  of  a  retardation  factor.  The 
basic  premise  for  using  retardation  factors  is  that  local  equilibrium  is 
a  valid  assumption.  Under  equilibrium  conditions,  a  retarded  solute 
will  move  through  an  aquifer  at  a  speed  equal  to  the  speed  with  which  a 
conservative  solute  moves,  divided  by  the  retardation  factor.  The 
retardation  factor  (R)  may  be  determined  by  any  of  the  following  equa¬ 


tions: 


p.'  (retarded) 

£ » i 

'  (conservative) 


(3-54a) 


ve^  (conservative) 
v  (retarded) 


(3-54b) 


using  either  spatial  or  temporal  definitions  for  vfi££,  and 

ujoo  (conservative) 


u[oo  (retarded) 


(3-54c) 


For  the  local  equilibrium  model,  all  of  the  above  equations  will 
yield  identical  values.  However,  if  the  assumption  of  local  equilibrium 
la  invalid  (physical  nonequilibrium,  for  example),  although  the  temporal 
moments  will  yield  the  same  retardation  factor,  independent  of  sampling 
location,  the  spatial  moments  will  yield  different  retardation  factors, 
depending  on  the  time  at  which  the  spatial  data  were  obtained.  The 
"nonequilibrium"  retardation  factors  obtained  from  the  spatial  data  will 
be  an  Increasing  function  of  time,  eventually  approaching  the  equilib¬ 
rium  retardation  factor  value  at  long  times. 

Another  consequence  of  the  "equilibrium”  behavior  of  the  temporal 
distribution  first  moment  will  now  be  discussed.  Van  Genuchten  and 
Wlerenga  (1976)  showed  that,  for  small  values  of  a,  the  first-order  rate 
constant,  a  temporal  breakthrough  response  could  be  simulated  using  a 
local  equilibrium  model  with 


and  for  large  values  of  a,  the  response  could  be  simulated  by  a  local 
equilibrium  model  with 


vef  f 


v 

1+6 


(3-56) 


Goltz  and  Roberts  (1986)  presented  similar  results,  and  further  defined 
small  and  large  values  of  a  relative  to  an  advective  rate  constant  (v/ £) 
such  that  if  a  «  v/£,  Eq.  3-55  would  apply,  and  if  a  »  v/£,  Eq.  3-56 
would  apply.  In  chemical  engineering  research,  the  ratio  of  a  mass 
transfer  rate  (a)  to  an  advective  rate  (v/ £)  is  defined  as  a  Stanton 
number  (Cussler,  1984),  where: 


St 


a 

7v7Ty 


(3-57) 


Thus,  for  St  «  1,  Eq.  3-55  would  apply,  and  for  St  »  1,  Eq.  3-56  would 
apply.  The  validity  of  (3-55)  and  (3-56)  is  illustrated  in  Figure  3.8a. 
For  the  simulation  at  a  low  mass  transfer  rate,  a  ■  0.0003  d  ^  and  v/ £  * 
0.10  d-1  (St  -  0.003).  Thus,  St  «  1,  and  the  local  equilibrium  model 
with  vg£f  ■  v  approximates  the  first-order  rate  response  quite  ade¬ 
quately.  For  the  simulation  at  a  high  mass  transfer  rate,  o  *  0.75  d 
and  v/ £  -  0.10  d”1  (St  -  7.5).  The  local  equilibrium  model  with 


-1 


'eff  “  l+B 


approximates  the  first-order  rate  response.  Interestingly,  an  examina¬ 
tion  of  the  temporal  moments  does  not  reveal  this  equivalence,  as  it  was 
found  that  the  first  temporal  moments  are  independent  of  the  mass  trans¬ 
fer  rate.  This  is  due  to  the  existence  of  a  very  small,  long  tail  in 
the  temporal  response,  which  causes  the  first  temporal  moment  to  remain 
constant  even  as  o  ♦  0.  Thus,  in  Figure  3.8a,  the  first-order  rate 
model  breakthrough  responses  for  the  high  and  low  mass  transfer  rates 
have  the  same  first  moment,  contrary  to  the  visual  impression  given  by 
the  figure.  Figure  3.8b,  which  depicts  the  tailing  of  the  high  and  low 
mass  transfer  rate  simulations,  indicates  how  the  long  tail  of  the  low 
rate  simulation  could  cause  the  two  simulations  shown  in  Figure  3.8a  to 
have  the  same  first  moment.  Although  mathematically  real,  the  tail 
which  leads  to  this  equality  may  be  below  detection  limits  (visual  and 
analytical),  and  it  is  reasonable  to  approximate  the  first-order  model 
response  at  small  a  using  the  local  equilibrium  model  with  vgff  as 
defined  in  Eq.  3-55.  The  equations  for  the  first  spatial  moment  provide 
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.  First-order  rate  model  simulations  for  high  (St  ■  7.5)  and 
low  (St  *  0.003)  values  of  the  first-order  mass  transfer 
rate  constant.  a)  Equilibrium  model  approximations,  and 
b)  Comparison  of  tailing  at  high  and  low  rates. 
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a  mathematical  justification  for  this  approximation.  Differentiate  the 
first  spatial  moment  expression  in  Table  3.6  to  find 


where 


_  dyioO  _  v(l  +  B2A)  ,  VBA(2  +  at  -  ate) 
eff  dt  TTTWTW  (1  +  BA)2 


(3-58) 


A  .  e-at(l+6) 


The  synoptic  sampling  time,  t,  in  Eq.  3-58  may  be  replaced  by  the  term 
l/v ,  to  obtain 


veff 


v(l  +  B2A) 

(1  +  B)(l  +  BA)  T 


vBA(2  + 

(1  +  BA)2 


(3-59) 


where 


A  .  e(-a£/v)(l+B) 

From  Eq.  3-59,  it  is  seen  that  indeed  for  a  «  v/l  (St  «  1) 

A  ♦  1  and  veff  ♦  v 

whereas  for  a  »  v/l  (St  »  1) 

A  +  0  <*nd  vef f  ♦  TTb 

These  approximations  could  not  be  obtained  by  an  examination  of  the 
first  temporal  moment. 

Immobile  region.  The  effective  velocity  of  the  immobile  solute 
distributions  based  on  temporal  moments  behaves  in  accordance  with 
Eq.  3-33.  At  small  distances  (£  «  v/[a(l+8)]  or  i  «  v/[a(l+8)]), 

vnef f  "  al  or  vnef f  “  <3-60> 

while  at  large  distances,  vneff  changes  (it  either  Increases  or  de¬ 
creases,  depending  on  parameter  values)  to 

vnef f  “  TTb  <3’61> 

As  discussed  above,  the  effective  velocity  of  the  immobile  region 
plume  calculated  from  the  spatial  moments  changes  from  the  small  time 
value  of 


for  the  first-order  rate  model,  and 


vnef f  *  HT  (3-63) 

for  the  diffusion  models,  to  the  large  time  value  of 

vnef f  "  ITb  (3-64) 


Second  moment 


The  value  of  Deff  calculated  from  the  temporal  moments  remains  con¬ 
stant  over  the  sampling  distance,  with  values  of 


Def f  “  Dx 


(3-65) 


for  the  local  equilibrium  model  and 

D 


x  v2B 

eff’1+6  o(l+B)3 


Def f  l  +  e  + 


v2  8 


(3-66) 


a(l  +  B)' 


for  the  physical  nonequilibrium  models. 

For  the  physical  nonequilibrium  models,  the  effective  dispersion 
coefficient  calculated  from  the  spatial  moments  increases  from 


Deffx  "  Dx 


(3-67) 


at  small  times,  to 


Def  fx 
Deffx 


3x  +  v2B 


1  +  B  '  .<1  ♦  6)3 


5x  +  V2B 


1  +  3  a(l  +  6)3 


(3-68) 


at  large  times.  Dejfx  for  the  local  equilibrium  models  remains  constant 
at  ®effx  “  ®x' 

Notice  that  all  the  large  time  parameter  values  obtained  using  the 
first  and  second  spatial  moments  are  equivalent  to  the  values  obtained 
using  the  temporal  moments. 


Table  3.9  summarizes  the  parameter  values  obtained  at  small  and 
large  values  of  sampling  distance  and  time,  based  on  the  temporal  and 
spatial  moments  respectively. 

SUMMARY 

One-  and  three-dimensional  forms  of  Arls'  method  of  moments  were 
developed  and  used  to  analyze  the  temporal  and  spatial  moment  behavior 
of  concentration  distributions  obtained  using  both  equilibrium  and 
physical  nonequilibrium  solute  transport  models. 

It  was  shown  that  although  the  zeroth  and  first  temporal  moments 
are  Independent  of  the  rate  of  mass  transfer  between  the  mobile  and 
immobile  regions,  all  the  spatial  moments  are  dependent  on  the  mass 
transfer  rate.  One  Implication  of  this  is  that  the  retardation  factors 
calculated  from  breakthrough  responses  should,  at  least  mathematically, 
be  Independent  of  the  distance  to  the  sampling  well  and  the  rate  of  mass 
transfer.  On  the  other  hand,  the  retardation  factors  calculated  from 
spatial  data  can  be  expected  to  increase  with  time  under  conditions  of 
physical  nonequilibrium,  with  the  rate  of  increase  a  function  of  the 
mass  transfer  rate. 

The  temporal  and  spatial  moment  behavior  of  the  solute  within  the 
immobile  region  was  also  examined.  It  was  demonstrated  that  the  zeroth 
spatial  moment  increases  from  zero,  indicating  no  solute  mass  within  the 
immobile  region  at  t  ■  0,  to  a  constant  value  at  long  time.  The  effec¬ 
tive  velocity  of  the  immobile  region  plume,  calculated  from  the  spatial 
moment,  changes  (it  may  increase  or  decrease)  from  a  fraction  of  the 
initial  mobile  plume  velocity  at  t  ■  0,  to  a  constant  value  at  long 
times,  with  the  large-time  value  equal  to  the  large-time  value  of  the 
mobile  plume  velocity. 

Finally,  it  was  shown  that  for  infinite  boundary  conditions,  the 
first  temporal  moment  simulated  using  a  one-dimensional  model  would  be 
greater  than  the  first  temporal  moment  using  a  three-dimensional  model 
with  the  same  velocity,  owing  to  the  effect  of  one-dimensional  versus 
three-dimensional  dispersion. 
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CHAPTER  4 
MODEL  COMPARISONS 


In  the  preceding  chapter,  the  spatial  and  temporal  moments  obtained 
using  different  transport  models  were  compared.  These  moment  formula¬ 
tions  will  now  be  used  to  define  equivalent  parameters  for  the  different 
transport  models,  and  to  examine  the  similarities  and  differences  of 
breakthrough  responses  (which  typically  would  be  measured  in  a  field 
situation)  simulated  using  these  equivalent  parameters.  In  contrast  to 
earlier  studies,  which  concentrated  on  equivalence  among  one-dimensional 
models,  this  work  will  deal  with  the  three-dimensional  models  discussed 
in  the  previous  chapters.  The  reason  for  this  analysis  of  model  equiva¬ 
lence  is  to  determine  if  differentiation  among  different  models,  based 
on  breakthrough  data,  is  feasible.  A  related  question,  which  will  also 
be  addressed  in  this  chapter,  concerns  the  differences  in  breakthrough 
responses  owing  to  different  assumptions  regarding  boundary  conditions. 


MODEL  EQUIVALENCE 

The  concept  of  equivalence  between  local  equilibrium  and  physical 
nonequilibrium  models  has  been  discussed  in  the  literature.  Baker 
(1977),  De  Smedt  and  Wierenga  (198A),  and  Valocchi  (1985a)  showed  that 
in  one  dimension,  for  Peffl  -  vi/Dx  >  1000,  the  first-order  rate  model 
could  be  approximated  by  a  local  equilibrium  model,  with  an  effective 
dispersion  coefficient  of: 


1  +  B 


a(l  +  6)' 


(4-1) 


Rao  et  al.  (1980),  and  Valocchi  (1985a)  showed  for  diffusion  into 
spheres,  and  Passloura  (1971)  showed  for  diffusion  into  cylinders  as 
well  as  into  spheres,  that  at  large  Pem  the  diffusion  models  could  be 
approximated  by  a  local  equilibrium  model  with  an  effective  dispersion 
coefficient  of: 


a(l  +  8)' 


(4-2) 


In  the  subsequent  discussion,  it  will  be  shown  that  Eq.  4-2  can  be  used 
to  approximate  diffusion  into  layers  also.  Valocchi  (1985a)  showed  that 


an  effective  velocity  for  a  local  equilibrium  model  equivalent  to  the 
physical  nonequilibrium  models  is: 

veff  *  rrr  (4_3) 

Equations  4-1  through  4-3  may  also  be  obtained  using  the  concept  of 
temporal  moment  equivalence.  From  Table  3.2  of  the  previous  chapter,  it 
is  seen  that  the  effective  parameter  values  expressed  in  Eqs.  4-1  through 
4-3  are  the  same  as  the  parameter  values  obtained  by  considering  the 
breakthrough  responses  of  the  one-dimensional  models  at  high  Peclet 
numbers,  and  by  setting  the  first  and  second  temporal  moments  of  the 
equilibrium  and  nonequilibrium  models  equal  to  each  other.  Inspecting 
the  first  moment  expressions  in  Table  3.2,  it  is  seen  that  at  high 
Peclet  numbers,  Eq.  4-3  can  be  used  in  equating  an  equilibrium  model 
effective  velocity  to  all  of  the  physical  nonequilibrium  models,  regard¬ 
less  of  whether  a  geometrical  diffusion  or  first-order  rate  model  is 
formulated. 

Comparing  the  second  moment  expressions  in  Table  3.2,  and  setting 
the  equilibrium  and  diffusion  model  expressions  equal  to  each  other,  at 
high  Peclet  numbers,  shows  Eq.  4-2  applicable  to  spherical,  cylindrical, 
and  layered  diffusion  models.  In  the  same  way,  Eq.  4-1  is  obtained  by 
setting  the  equilibrium  model  expression  for  the  second  temporal  moment 
equal  to  the  first-order  rate  model  expression  for  the  same  moment.  To 
define  equivalence  between  first-order  rate  and  diffusion  model  parame¬ 
ters,  it  is  necessary  merely  to  set  the  second  temporal  moment  formula¬ 
tions  of  the  two  models,  shown  in  Table  3.2,  equal  to  each  other,  to 
find: 

o  *  a  (4-4) 

Note  that  the  validity  of  Eq.  4-4  does  not  depend  on  the  Peclet  number. 

Parker  and  Valocchi  (1986)  used  the  moment  analysis  approach  de¬ 
scribed  above  to  define  equivalent  equilibrium,  first-order  rate,  and 
spherical  diffusion  model  parameters.  Here,  their  analysis  is  extended 
to  include  cylindrical  and  layered  diffusion  models. 

Van  Genuchten  (1985)  used  an  empirical  method  to  obtain  shape 
factors  which  could  be  used  to  define  equivalent  rate  constants  for  the 
various  physical  nonequilibrium  models.  Similar  factors  can  be  obtained 
using  Eq.  4-4.  Table  4.1  compares  equivalent  rate  constants  obtained  using 


TABLE  4.1 

EQUIVALENT  RATE  CONSTANTS  FOR  NONEQUILIBRIUM  MODELS 


Diffusion  Models 

First-Order  - 

Rate  Model  Spherical  Cylindrical  Layered 


Empirical 


* 
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Moment  Analysis 
(Eq.  4-4) 


*Van  Genuchten  (1985). 
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van  Genuchten' s  (1985)  empirical  formulae  with  rate  constants  derived 
using  the  moment  analysis  described  above. 

Van  Genuchten  (1985)  and  Parker  and  Valocchl  (1986)  provide  compar¬ 
isons  of  the  one-dimensional  breakthrough  responses  of  the  equilibrium 
and  nonequilibrium  models.  In  this  work,  it  is  shown  that  their  find¬ 
ings  are  also  applicable  in  three  dimensions. 

With  regard  to  the  physical  nonequilibrium  models,  Parker  and 
Valocchi  (1986)  considered  only  the  first-order  rate  and  spherical  dif¬ 
fusion  models.  However,  as  was  shown  by  van  Genuchten  (1985)  and  is 
apparent  from  Figure  3.4c,  these  two  model  types  bracket  the  range  of 
behavior  of  the  other  geometrical  diffusion  models.  Therefore,  the  fol¬ 
lowing  three-dimensional  analysis  will  follow  Parker  and  Valocchl  (1986) 
and  be  limited  to  discussing  the  first-order  rate  and  spherical  diffu¬ 
sion  models.  Parker  and  Valocchi  (1986)  also  found  that  the  deviations 
between  models  are  greater  for  an  instantaneous  source  than  for  continu¬ 
ous  or  pulse  injections.  Hence,  the  following  discussion  will  focus  on 
the  breakthrough  responses  to  an  instantaneous  point  source. 

The  instantaneous  point  source  analysis  which  will  be  discussed  in 
this  work  is  applicable  to  sampling  wells  far  from  a  finite  source,  in 
the  line  of  advective  transport.  Parameter  values  used  for  the  analysis 
are  listed  in  Table  4.2.  The  parameter  y  is  a  measure  of  the  ratio  of 
diffusive  to  advective  rate  constants.  For  the  first-order  rate  model: 
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TABLE  4.2 

PARAMETER  VALUES  USED  IN  FIGURES  4.1  THROUGH  4.3 


Y 

6 

R 

Peef  f 

Pem 

Figure  4.1 

30.0 

1 

1 

27.3 

50 

0.3 

1 

1 

0.593 

50 

Figure  4.2 

0.3 

1 

1 

0.593 

50 

0.3 

0.0526 

1 

18.3 

50 

Figure  4.3 

0.3 

1 

1 

0.593 

50 

0.3 

1 

50 

18.7 

50 

and  for  the  diffusion  models: 

» ■  wzy  <4-6> 

Note  that  in  the  case  of  the  first-order  rate  model  (Eq.  4-5),  y  is  equiv¬ 
alent  to  the  Stanton  number  defined  in  Chapter  3.  In  the  case  of  the 
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ire  4.2.  Effect  of  6  on  simulated  breakthrough  responses  of  the 
equilibrium,  first-order  rate,  and  spherical  diffusion 
models. 
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Figure  4.3.  Effect  of  R  on  simulated  breakthrough  responses  of  the 
equilibrium,  first-order  rate,  and  spherical  diffusion 
models. 


modulus  (Crittenden  et  al.,  1986).  Nonequilibrium  model  parameters  were 
chosen  to  provide  a  range  of  behavior.  The  equivalent  equilibrium  model 
parameter  values  were  calculated  using  Eqs.  4-1  through  4-3.  The  equiv¬ 
alent  equilibrium  model  parameter  values  are  contained  in  the  effective 
Peclet  number,  Pgff  -  veffi/Deff. 

Both  van  Genuchten  (1985)  and  Parker  and  Valocchi  (1986)  found  that 
for  large  values  of  y,  the  ratio  of  the  diffusive  rate  constant  to  the 
advective  rate  constant  (i.e.,  approach  to  equilibrium),  the  simulations 
of  the  local  equilibrium,  first-order  rate,  and  diffusion  models  con¬ 
verged  upon  each  other.  Figure  4.1,  which  plots  normalized  concentration 

A 

versus  normalized  time  (T),  where 


demonstrates  this  for  a  three-dimensional  simulation.  For  high  values 
of  y,  the  equilibrium  and  nonequilibrium  models  produce  similar  break¬ 
through  reponses,  whereas  for  low  values  of  y»  the  equilibrium  and  non- 
equilibrium  responses  diverge. 

Figure  4.2  shows  the  effect  of  B,  the  solute  capacity  ratio  of 
Immobile  to  mobile  regions,  upon  simulated  breakthrough  responses.  As  B 
approaches  zero  (i.e.,  where  very  little  solute  is  associated  with  the 
immobile  region),  the  nonequilibrium  responses  approach  the  equilibrium 
responses.  For  higher  values  of  B>  which  implies  more  solute  associated 
with  the  Immobile  region,  the  deviation  from  equilibrium  Increases. 
Parker  and  Valocchi  (1986)  showed  a  similar  effect  in  one  dimension. 

The  effect  of  increasing  the  retardation  factor,  R,  upon  the  simu¬ 
lated  responses  is  shown  in  Figure  4.3.  As  can  be  seen,  increasing  R 
decreases  the  deviation  from  equilibrium. 

Figures  4.1  through  4.3  show  that  for  certain  parameter  values, 
even  if  Pe,,  <  1000,  Eqs.  4-1  through  4-3  are  applicable.  This  was  also 
shown  by  Parker  and  Valocchi  (1986).  The  figures  also  indicate  that  for 
other  parameter  values,  the  differences  between  the  geometrical  diffu¬ 
sion  models  and  the  equivalent  first-order  rate  model  may  be  significant. 
This  observation  was  also  made,  for  one-dimensional  formulations,  by  van 
Genuchten  (1985)  and  Parker  and  Valocchi  (1986).  Van  Genuchten  (1985) 
noted,  however,  that  considering  the  uncertainty  in  the  many  parameters 
needed  in  the  two-region  models,  the  first-order  rate  model  acceptably 


approximated  the  responses  of  the  more  complex  diffusion  models.  In  the 
above  analysis,  the  differences  in  model  responses  were  maximized  by 
choosing  an  instantaneous  point  source  and  comparing  the  two  most  dis¬ 
similar  models.  Nevertheless,  the  differences  between  the  responses 
were  not  very  great.  This  observation  supports  van  Genuchten's  (1985) 
conclusion  that,  owing  to  parameter  uncertainty,  attempts  to  differen¬ 
tiate  between  the  two  types  of  models  based  on  field  data  would  be 
difficult. 

Figure  4.4  offers  further  evidence  in  support  of  this  conclusion. 
As  Figures  4.1  through  4.3  showed,  for  certain  parameter  values,  the 
differences  between  the  first-order  rate  and  spherical  diffusion  models 
could  be  significant.  However,  with  relatively  minor  changes  in  the 
first-order  rate  model  parameter  values,  the  responses  of  the  two  models 
can  be  made  to  be  almost  identical.  In  Figure  4.4,  v  and  Dx  were 
changed  by  10%  from  the  Figure  4.1  through  4.3  values,  and  the  first- 
order  rate  constant,  a,  was  calculated  using  van  Genuchten's  (1985) 
empirical  formula  (see  Table  4.1)  instead  of  using  Eq.  4-4.  As  the 
figure  shows,  the  responses  of  the  two  models  are  quite  similar. 
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Figure  4.4.  Comparison  of  simulated  breakthrough  responses  of  the 
first-order  rate  and  spherical  diffusion  models. 


Parker  and  Valocchi  (1986)  showed  that  for  certain  parameter 
values,  it  is  possible  for  a  local  equilibrium  model  simulation  to 
deviate  even  less  than  a  first-order  rate  model  from  a  diffusion  model 
simulated  breakthrough  response.  Thus,  the  likelihood  of  successfully 
differentiating  among  the  three  types  of  models  based  on  field  break¬ 
through  response  data  appears  to  be  small. 

Recall,  however,  from  the  discussion  in  the  previous  chapter  that 
there  is  a  qualitative  difference  between  the  spatial  moment  predictions 
of  the  equilibrium  and  nonequilibrium  models.  Thus,  from  spatial  data, 
it  may  be  possible  to  differentiate  between  equilibrium  and  nonequilib¬ 
rium  models,  although  as  can  be  seen  from  Figure  3.4,  differentiating 
between  the  different  nonequilibrium  model  formulations  remains  diffi- 


INFINITE  AND  SEMI-INFINITE  BOUNDARY  CONDITIONS 

To  examine  the  effect  of  infinite  versus  semi-infinite  boundary 
conditions  on  the  breakthrough  responses,  one-dimensional  breakthrough 
simulations  of  semi-infinite  and  infinite  versions  of  the  first-order 
rate  and  spherical  diffusion  models  will  be  compared  in  this  section. 

Lindstrom  and  Narasimhan  (1973)  derived  the  first-order  rate  model 
solution  for  an  initially  distributed  solute  slug  in  a  one-dimensional 
semi-infinite  medium.  The  analogous  solution,  for  the  following  initial 
and  boundary  conditions  corresponding  to  an  infinitesimal  point  source: 
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and  H(t,T)  is  as  defined  In  Eq.  2-8. 

Using  Appendix  A,  the  one-dimenslonal  version  of  Eq.  2-8  may  be 
written 


C^x.t)  ■  exp[-  -|-|-]G(x,t)  +  /  H(t,x )G(x,x )dt 
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(4-10) 
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Equation  4-10  represents  the  response  to  an  Instantaneous  point  source 
In  a  one-dimensional  infinite  medium. 

Wakao  and  Kaguei  (1982)  examined  the  differences  in  using  either 
semi-inf inite  or  infinite  boundary  conditions  in  a  model  which  included 
both  advective/disperslve  transport  and  diffusion  into  zones  of  immobile 
water.  They  found  that  deviations  between  the  two  solutions  were  depen¬ 
dent  on  the  dimensionless  group,  Pem  -  ^m^mx)’  Fi8ure  4.5  shows 
dimensionless  concentration  versus  time  solutions  of  the  first-order 
rate  model  in  both  a  semi-infinite  (4-9)  and  an  infinite  (4-10)  medium, 
for  various  values  of  Pem.  For  low  values  of  Pem,  differences  between 
the  two  solutions  become  significant.  At  low  values  of  Pem,  the  disper¬ 
sive  term  in  the  advection/dispersion  equation  is  large  enough  to  cause 
differences  between  the  semi-infinite  solution,  which  does  not  allow  for 
dispersive  solute  flux  into  the  region  x  <  0,  and  the  infinite  solution, 
which  does  permit  upgradient  dispersion.  As  expected,  Figure  4.5  shows 
that  at  low  values  of  Pem,  the  infinite  solution  exhibits  more  spreading 
than  does  the  semi-infinite  solution.  At  high  values  of  Pem,  however 
(i.e.,  Pea  >  50),  the  differences  between  the  semi-infinite  and  infinite 
solutions  are  insignificant. 
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Figure  4.5. 


Comparison  of  semi-infinite  and  infinite  solutions  to  the 
first-order  rate  model:  <|>  *  0.90,  a  *  0.004  d~*,  v/ J l  * 
0.02  d_1,  and  R^  *  Rim  =  1.0. 


The  solution  of  the  spherical  diffusion  model  in  a  one-dimensional 
semi-infinite  medium  is  found  by  a  straightforward  application  of  the 
methods  of  Rasmuson  and  Neretnieks  (1980)  and  van  Genuchten  et  al. 
(1934).  The  solution,  for  an  Instantaneous  point  source,  is: 


Using  the  methods  of  Appendix  B,  the  one-dimensional  version  of 
Eq.  2-12  may  be  written: 


where,  again,  Zp  and  Zm  are  defined  in  Eq.  2-12.  Equation  4-12 
represents  the  response  to  an  Instantaneous  point  source  in  a  one¬ 
dimensional  infinite  medium. 


Figure  4.6  compares  dimensionless  concentration  versus  time 
solutions  of  the  spherical  diffusion  model  in  both  a  semi-infinite 
(4-11)  and  an  infinite  (4-12)  medium,  for  various  values  of  Pem.  The 
results  are  quite  similar  to  the  results  obtained  using  the  first-order 
rate  model.  At  high  values  of  Pem,  the  difference  between  the  semi¬ 
infinite  and  infinite  solutions  decreases. 

CONCLUSIONS 

In  this  chapter,  the  equivalence  between  the  local  equilibrium, 
first-order  rate,  and  geometrical  diffusion  models  was  evaluated.  It  was 
shown  that  the  equivalence  relationships  derived  using  one-dimensional 
model  formulations  of  the  various  models  are  also  applicable  in  three 
dimensions,  and  that  breakthrough  responses  simulated  using  the  one-  and 
three-dimensional  models  show  a  similar  dependence  on  input  parameter 
values . 

Use  of  field  breakthrough  (temporal)  data  to  differentiate  between 
the  three  types  of  models  does  not  seem  feasible,  though  spatial  moment 
data  may  be  useful  in  differentiating  between  the  equilibrium  and  non¬ 
equilibrium  model  formulations. 
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Figure  4.6.  Comparison  of  semi-infinite  and  infinite  solutions  to  the 
spherical  diffusion  model:  4>  *  0.90,  De/b^  *  0.004  d“  , 
v/£  -  0.02  d  ,  and  *  Rim  «  1.0. 


It  was  also  demonstrated  that  the  solutions  for  semi-infinite  and 
infinite  boundary  conditions  are  similar  at  large  values  of  the  Peclet 
number.  This  is  convenient,  since  the  infinite  solution  is  simpler  and 
obtained  in  a  more  straightforward  manner  than  the  semi-infinite  solu¬ 
tion. 
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CHAPTER  5 

APPLICATION  TO  THE  INTERPRETATION  OF  DATA  FROM  A 
LARGE-SCALE  TRANSPORT  EXPERIMENT  UNDER  NATURAL  CONDITIONS 

The  discussion  thus  far  has  concentrated  on  the  development  and 
analysis  of  three-dimensional  physical  nonequilibrium  models.  In  this 
chapter,  spatial  and  temporal  data  obtained  from  a  large-scale  transport 
experiment  will  be  interpreted  using  such  models.  The  physical  nonequi¬ 
librium  models,  with  parameters  obtained  from  direct  measurement,  liter¬ 
ature  correlations,  and  laboratory  experiments,  will  be  used  in  an 
attempt  to  simulate  the  spatial  and  temporal  behavior  of  the  solute 
distributions  which  were  observed  over  the  course  of  the  experiment.  In 
addition,  the  complementarity  of  the  results  from  spatial  and  temporal 
sampling  will  be  assessed  in  two  ways.  First,  a  direct  comparison  will 
be  made  of  spatial  and  temporal  results  obtained  at  comparable  time/dis¬ 
tance  scales.  Second,  equilibrium  and  physical  nonequilibrium  model 
parameters  obtained  from  spatial  sampling  will  be  used  to  simulate 
temporal  response  data.  Finally,  a  brief  review  will  be  made  of  various 
alternative  models  which  may  be  used  in  interpreting  the  field  results. 
This  review  will  include  an  examination  of  the  spatial  moment  behavior 
predicted  by  some  of  the  alternative  models,  a  subject  that  previously 
has  not  been  evaluated  systematically. 

PROJECT  BACKGROUND 

A  large-scale,  long-term  field  experiment  to  study  natural  gradient 
transport  of  solutes  in  groundwater  was  conducted  by  a  group  of  investi¬ 
gators  from  Stanford  University  and  the  University  of  Waterloo,  at  a 
site  in  Borden,  Ontario  (Mackay  et  al.,  1986;  Freyberg,  1986;  Roberts  et 
al.,  1986).  Well-defined  initial  conditions  were  achieved  by  Injecting 
approximately  12  nr  of  a  solution  containing  known  masses  of  two  inor¬ 
ganic  tracers  (chloride  and  bromide)  and  five  halogenated  organic 
compounds  (bromoform,  carbon  tetrachloride,  tetrachloroethylene,  1,2- 
dichlorobenzene,  and  hexachloroe thane ) .  The  transport  of  the  organic 
solutes  was  monitored  over  a  two-year  period  with  a  dense,  three- 
dimensional  array  of  more  than  4000  sampling  points.  Figure  5.1  depicts 
the  sampling  network. 
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Figure  5.1.  Locations  of  multilevel  sampling  and  injection  wells  as  of 
January  1986:  a)  Plan  view,  and  b)  Approximate  vertical 
distribution  of  sampling  points  (+)  projected  onto  cross 
section  AA*  (vertical  exaggeration  -  4.6). 


Two  types  of  sampling  were  conducted.  Spatial  data  were  obtained 
from  synoptic  or  "snapshot"  sampling  sessions ,  each  of  which  measured 
the  three-dimensional  spatial  distribution  of  solute  concentration  at  a 
particular  point  in  time.  Temporal  data  were  obtained  by  measuring  the 
solute  concentration  at  a  relatively  high  sampling  frequency  (i.e., 
daily  in  the  early  stages,  biweekly  or  weekly  later  on)  at  a  few  pre¬ 
selected  sampling  points.  Tables  5.1  and  5.2  summarize  the  synoptic  and 
time  series  monitoring  programs.  Details  of  the  experimental  design  and 
implementation  may  be  found  in  Mackay  et  al.  (1986). 

The  aquifer  in  which  the  experiment  was  conducted  is  unconfined, 
consisting  of  horizontally  bedded  fine-  to  medium-grained  sand  (Sudicky 
et  al.,  1983).  The  hydrogeology  and  geochemistry  of  the  study  area  have 
been  described  by  MacFarlane  et  al.  (1983)  and  Dance  (1980). 


TABLE  5.1 


SUMMARY  OF  SYNOPTIC  MONITORING  PROGRAM 


Date 

Days  Since 
Injection 

Solutes 

Tracers  Organics 

Number  of  Samples 
Analyzed 

08/24/82 

1 

X 

X 

392 

09/01/82 

9 

X 

X 

419 

09/08/82 

16 

X 

X 

408 

09/21-22/82 

29 

X 

X 

629 

10/05-06/82 

43 

X 

X 

671 

10/25-26/82 

63 

X 

X 

700 

11/16-17/82 

85 

X 

X 

712 

05/09-11/83 

259 

X 

1219 

06/22/83 

303 

X 

233 

07/19-20/83 

330 

X 

1150 

07/21/83 

332 

X 

362 

09/07-08/83 

380 

X 

839 

09/08-09/83 

381 

X 

496 

10/04/83 

407 

X 

949 

10/26-28/83 

429 

X 

1883 

11/28/83 

462 

X 

1343 

05/17/84 

633 

X 

1122 

05/31-06/02/84 

647 

X 

958 

08/01-02/84 

709 

X 

1119 

06/26-28/85 

1038 

X 

1205 

TABLE  5.2 

SUMMARY  OF  TIME  SERIES  MONITORING  PROGRAM 


Sample  Point  (x,y,z) 

(m) 


Duration 


Number  of  Samples 
Collected  as  of  1/1/86 


Near- 

Field 


Mid- 

Field 


Far- 

Field 


2.5, 

0.0,  - 

-3.20 

2.5, 

1.25, 

-3.62 

5.0, 

0.0,  - 

-3.26 

10.0, 

4.6, 

-3.88 

10.0, 

4.6, 

-4.48 

13.1, 

4.05, 

,  -3.42 

13.1, 

4.05 

,  -3.72 

18.0, 

9.0, 

-4.13 

18.0, 

9.0, 

-4.73 

21.0, 

9.0, 

-4.17 

21.0, 

9.0, 

-4.77 

24.0, 

9.0, 

-4.76 

Aug.  82  -  Dec. 
Aug.  82  -  Dec. 
Aug.  82  -  Dec. 

Nov.  83  -  Jun. 
Nov.  83  -  Jun. 
Jul.  84  -  Jun. 
Jul.  84  -  Jun. 

Mar.  83  -  Jun. 
Mar.  83  -  Jun. 
Mar.  83  -  Jun. 
Mar.  83  -  Jun. 
Mar.  83  -  Nov. 


83 

188 

83 

188 

83 

183 

85 

45 

85 

27 

85 

32 

85 

31 

85 

119 

85 

121 

85 

117 

85 

117 

83 

78 

Time  Series  Data 

Concentration  responses  at  each  of  the  twelve  high  frequency  sampl¬ 
ing  points  are  shown  in  Appendix  E.  Chloride  data  are  not  depicted,  in 
the  interest  of  clarity,  since  the  chloride  and  bromide  observations  are 
nearly  Indistinguishable  (Freyberg,  1986).  For  the  purpose  of  this 
study,  1,2-dichlorobenzene  and  hexachloroethane  data  have  also  been 
omitted.  As  discussed  in  Roberts  et  al.  (1986),  both  of  these  compounds 
behaved  anomalously.  Hexachloroethane  concentrations  declined  to  non- 
quantlfiable  levels  by  the  end  of  the  first  three  months  of  the  experi¬ 
ment,  and  the  concentration  of  1,2-dichlorobenzene  declined  precipitously 
at  one  of  the  three  near-field  sampling  wells  and  was  never  found  in 
significant  levels  thereafter.  Significant  concentration  levels  of 
1 ,2-dichlorobenzene  were  not  seen  at  either  of  the  other  two  near-field 
wells.  A  full  discussion  of  the  near-field,  early-time  behavior  of 
these  two  compounds  may  be  found  in  Roberts  et  al.  (1986).  This  study 
will  focus  on  those  compounds  that  remained  at  significant  concentration 
levels  after  the  first  three  months  of  monitoring. 

Appendix  E  also  shows  that  of  the  four  solutes  being  considered, 
every  solute  is  not  seen  at  every  sampling  point.  The  four  mid-field 
wells  were  sampled  specifically  for  tetrachloroethylene .  Bromide  and 
the  two  faster  moving  organics,  carbon  tetrachloride  and  bromoform,  had 
largely  passed  by  these  wells  before  sampling  commenced.  At  the  five 
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far-field  wells,  tetrachloroethylene  had  yet  to  arrive  upon  conclusion 
of  the  experiment.  At  the  most  distant  well,  only  bromide  was  seen, 
since  sampling  at  this  well  was  terminated  prior  to  the  arrival  of  the 
organic  solutes. 

Synoptic  Data 

Freyberg  (1986)  described  how  the  synoptic  concentration  data  were 
used  to  obtain  estimates  of  the  spatial  moments  of  the  solute  plume  dis¬ 
tributions.  Freyberg  (1986)  and  Roberts  et  al.  (1986)  present  estimates 
of  the  zeroth  and  first  moments  for  the  inorganic  and  organic  plumes, 
respectively.  Freyberg  (1986)  also  presents  estimates  of  the  second 
moment  for  the  inorganic  plumes.  Following  his  methods,  the  organic 
plumes'  second  moment  estimates  were  calculated  as  well  (Freyberg, 
1985).  In  this  work,  estimates  of  the  second  moment  are  presented  as 
(Jxx  and  Oyy,  the  principal  values  of  the  spatial  covariance  tensor. 
Freyberg  (1986)  noted  that  the  vertical  components  of  the  spatial 
covariance  tensor  could  not  be  distinguished  from  sampling  noise,  and 
that  the  vertical  thickness  of  the  tracer  plumes  remained  essentially 
constant  over  the  course  of  the  experiment.  Therefore,  the  principal 
vertical  component  of  the  covariance  tensor,  ozz,  has  not  been  consid¬ 
ered  and  will  be  assumed  to  be  negligible. 


Comparison  of  Time  Series  and  Synoptic  Results 


Roberts  et  al.  (1986)  demonstrated  that  the  results  of  time  series 
sampling  of  the  near-field  wells  (x  <_  5  m)  were  qualitatively  consistent 
with  synoptic  sampling  results.  In  this  section,  far-field  (x  >  15  m) 
time  series  sampling  results  for  bromoform  and  carbon  tetrachloride  are 
presented,  and  compared  with  synoptic  sampling  results,  to  further 
demonstrate  the  complementarity  of  spatial  and  temporal  transport  behav¬ 
ior.  The  comparison  is  made  in  terms  of  the  retardation  factor,  which 
is  calculated  using  temporal  or  spatial  first  moment  estimates,  and  the 
relative  mass  ratio,  which  is  obtained  from  temporal  or  spatial  zeroth 
moment  estimates.  The  following  definitions  are  used: 

1.  Retardation  factors  are  defined  by  Eqs.  3-55a  and  3-55c, 
based  on  temporal  and  spatial  first  moment  estimates, 
respectively. 
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2.  Relative  mass  is  defined  as  the  total  solute  mass  calcu¬ 
lated  from  the  temporal  or  spatial  zeroth  moment  estimate, 
normalized  by  the  total  solute  mass  initially  injected 
into  the  aquifer.  The  relative  mass  ratio  is  the  rela¬ 
tive  mass  of  a  retarded  solute  divided  by  the  relative 
mass  of  the  conservative  tracer.  The  relative  mass  cal¬ 
culated  from  the  temporal  zeroth  moment  is  proportional 
to  the  area  under  the  breakthrough  response  curve.  The 
relative  mass  calculated  from  the  spatial  zerotii  moment 
is  proportional  to  the  product  of  the  zeroth  moment  and 
the  retardation  factor.  As  noted  above,  for  this  analy¬ 
sis,  the  retardation  factor  at  a  particular  synoptic 
sampling  time  is  obtained  from  Eq.  3-55c.  The  reason  the 
mass  calculated  from  the  spatial  moment  is  proportional 
to  the  retardation  factor  while  the  mass  calculated  from 
the  temporal  moment  is  not,  goes  back  to  the  discussion 
in  Chapter  3.  When  determining  the  spatial  moment  at  an 
Instant  in  time,  only  the  aqueous  solute  concentration  is 
measured.  The  sorbed  solute  must  be  accounted  for,  in 
this  instance,  by  use  of  the  retardation  factor.  On  the 
other  hand,  when  determining  the  temporal  zeroth  moment, 
all  the  solute  passing  by  a  particular  sampling  point  is 
accounted  for  in  the  breakthrough  response  curve. 

In  order  to  graphically  compare  spatial  and  temporal  behavior,  it 
is  useful  to  plot  the  two  types  of  results  on  the  same  scale.  Spatial 
data  are  obtained  at  particular  sampling  times.  The  sampling  distances 
at  which  temporal  data  are  obtained  may  be  converted  to  equivalent  time 
values.  For  this  analysis,  the  first  moment  of  the  breakthrough  re¬ 
sponse  curve  of  the  solute  of  interest  at  a  sampling  point  is  defined  as 
the  equivalent  time  for  that  point. 

Figures  5.2  and  5.3  depict  the  relative  mass  ratios  and  the  retar¬ 
dation  factors  for  bromoform  and  carbon  tetrachloride  calculated  from 
temporal  and  spatial  moment  estimates.  For  bromoform,  the  decrease  in 
mass  and  the  increase  in  the  retardation  factor  which  was  found  by  an 
examination  of  the  spatial  data  (Roberts  et  al.,  1986)  is  also  clearly 
seen  by  comparing  near-  and  far-field  temporal  results.  The  obvious 
disappearance  of  bromoform  (Figure  5.2)  was  attributed  to  transformation 
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□  Estimated  from  spatial  data 


4-  Estimated  from  temporal  data 


□  Estimated  from  spatial  data 


4-  Estimated  from  temporal  data 


TIME  (doys) 


Figure  5 .2.  Comparison  of  the  bromoform  plume  behavior  determined  from 
spatial  and  temporal  data:  a)  Relative  mass  ratio,  and 
b)  Retardation  factor  versus  time. 
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Figure  5.3.  Comparison  of  the  carbon  tetrachloride  plume  behavior 
determined  from  spatial  and  temporal  data:  a)  Relative 

mass  ratio,  and  b)  Retardation  factor  versus  time. 
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by  Roberts  et  al.  (1986).  More  variability  is  evident  in  the  calcula¬ 
tions  of  the  relative  mass  for  carbon  tetrachloride,  though  again,  the 
increase  in  the  retardation  factor  over  the  course  of  the  experiment  is 
evident.  The  variability  of  the  relative  mass  ratios  obtained  from  the 
temporal  data  may  reflect  local  variations  in  aquifer  properties  along 
the  flow  paths  to  the  Individual  sampling  points  (Roberts  et  al.,  1986). 
Nevertheless,  there  is  good  general  agreement  between  the  results  of  the 
spatial  and  temporal  analyses. 


MODEL  APPLICATION 

One  goal  of  modeling  is  to  predict  real  behavior  using  indepen¬ 
dently  obtained  parameter  values.  In  this  section,  an  attempt  is  made 
to  predict  zeroth  and  first  spatial  moment  behavior  using  first-order 
rate  model  parameters  equivalent  to  the  spherical  and  layered  diffusion 
model  parameters.  Model  parameter  values  are  estimated  based  on  direct 
measurement,  literature  correlations,  and  soil  adsorption  isotherms 
measured  in  the  laboratory. 


Spherical  Diffusion  Model 


Based  on  the  approximately  spherical  shape  of  the  aquifer  material, 
it  is  reasonable  to  assume  solute  diffusion  into  spherical  grains  may  be 
influencing  transport.  Table  5.3  lists  parameter  values  for  use  in  a 
spherical  diffusion  model.  The  parameters  in  Table  5.3  were  obtained 
Independently  of  the  experimental  data.  The  mass  of  each  solute  ini¬ 
tially  Injected  into  the  aquifer  (M^)  was  known  (Mackay  et  al.,  1986). 
Values  for  bulk  density  (p),  and  total  porosity  (6)  were  obtained  by 
direct  measurement  (Mackay  et  al.,  1986)  as  was  the  range  of  median 
grain  sizes  (b)  (O'Hannesln,  1981).  The  value  for  the  pore  water 
velocity  (v0)  was  based  on  estimates  of  hydraulic  conductivity  and 
hydraulic  gradient  (Mackay  et  si.,  1986).  The  ratio  of  mobile  to  total 
water  content  ($)  was  calculated  using  the  following  expression: 


(5-1) 


where  ea  is  the  intragranular  porosity.  Intragranular  porosity  was 
measured  by  mercury  porosimetry  (Demond,  1984).  The  distribution  coeffi¬ 
cients  (Kj)  for  the  various  organic  compounds  and  Borden  aquifer  material 


TABLE  5.3 

PARAMETER  VALUES  FOR  USE  IN  A  SPHERICAL  DIFFUSION  MODEL 


Bulk  soil  density  (p) 

1.81  g/cm3  i 

Total  porosity  (0) 

0.33 

4 

Intraparticle  porosity  (e_) 

0.01 

•  I 

Ratio  of  mobile  to  total  water  content 

(4>) 

0.979  i 

Pore  water  velocity  (v  ) 

0.079  m/d  J 

Fraction  of  "mobile"  sorption  sites  (f) 

0.0 

1 

Median  grain  radius  (b) 

0.035-0.35  mm  j 

Carbon 

Tetra-  1 

Bromo- 

Tetra- 

chloro-  J 

Bromide 

form 

chloride 

ethylene  1 

Distribution  coeff. 

} 

4 

(Kd),  cm3/ g 

0.0 

0.17 

0.17 

0.48  i 

Mass  of  Injected 
solute  (Mp ,  g 

Liquid  diffusion 

3,870 

0.38 

0.37 

0.36  i 

coeff.  (DQ),  m2/d 

1.74X10-4 

7.17xl0-5 

7.34xl0-5 

6.83xl0~5  J 

Effective  intraparticle 

1 

• 

diffusion  coeff. 

(Dp,  m2/d 

1.74xl0“6 

7.17xlO-7 

7.34xl0-7 

6 .83  xlO-7  ■ 

Equivalent  first-order 

1 

rate  constant  (a'), 
d-1 

1.44-144 

0.594-59.4 

0.608-60.8 

0.566-56.6  S 

were  measured  by  Curtis  et  al.  (1986)  In  a  series  of  batch  sorption  ex¬ 
periments.  Bromide,  the  Inorganic  tracer,  is  assumed  to  be  nonsorbing. 
The  liquid  diffusion  coefficients  (DQ)  for  the  organic  compounds  were 
estimated  using  the  Wilke  and  Chang  (1955)  correlation.  The  liquid 
diffusion  coefficient  for  bromide  was  calculated  using  a  method  de¬ 
scribed  by  Cussler  (1984)  for  estimating  the  diffusion  coefficient  of 
electrolytes.  Since  Intraparticle  diffusion  is  likely  to  be  hindered  by 
the  Interaction  of  the  solute  and  the  sorbent  grain,  the  effective 
Intraparticle  diffusion  coefficient  (Dg)  is  some  fraction  of  the  liquid 
diffusion  coefficient.  Defining  tortuosity  (x)  such  that: 


Given  the  measured  value  of  e&  of  0.01,  x  “  100.  Equation  5-2  was  used 
to  calculate  the  effective  intraparticle  diffusion  coefficients  for  the 
compounds  being  considered.  Finally,  a  value  for  f,  the  fraction  of 
sorption  sites  adjacent  to  regions  of  mobile  water  is  required  for  model 
input.  Nkedi-Kizza  et  al.  (1982)  assumed  that  for  an  aggregated  soil, 
virtually  all  the  sorption  occurs  within  the  aggregate,  so  that  there 
are  no  sorption  sites  in  direct  contact  with  mobile  water.  Following 
Nkedi-Kizza  et  al.  (1982),  it  will  be  assumed  that  sorption  occurs 
either  within  the  aquifer  material  grains  or  within  aggregates,  so  that 
f  ■  0.0  for  this  model. 

It  is  now  possible  to  apply  the  moment  analysis  formulae  of  Table 
4.1  to  obtain  first-order  rate  constants  (o')  equivalent  to  the  spheri¬ 
cal  diffusion  rates  which  can  be  calculated  using  the  values  of  Dg,  b, 
$,  and  6  in  Table  5.3.  This  conversion  to  a  first-order  rate  model  is 
chosen  because  of  the  convenience  and  computational  efficiency  of  the 
first-order  rate  model  as  compared  to  the  diffusion  models.  In  Chap¬ 
ter  4,  it  was  shown  that  simulations  obtained  from  each  of  the  two  model 
types  are  quite  similar,  and  probably  Indistinguishable  considering  the 
accuracy  of  the  experimental  data  and  the  accuracy  with  which  the  input 
parameters  are  estimated.  The  range  in  values  for  a'  is  due  to  the 
range  in  median  grain  sizes. 

Before  using  the  parameter  values  estimated  above  in  the  first- 
order  rate  model,  and  comparing  simulations  with  experimental  data,  it 
is  useful  to  discuss  the  relevance  of  the  boundary  and  initial  condi¬ 
tions  assumed  in  the  model  to  the  experimental  conditions.  The  assump¬ 
tion  of  Infinite  boundary  conditions  would  seem  applicable  in  describing 
the  experimental  conditions  at  Borden,  in  which  the  medium  is  unbounded 
upgradient  from  the  point  of  introduction,  as  well  as  downgradient  from 
the  observation  points.  In  any  event,  as  was  discussed  in  Chapter  4, 
except  for  observation  points  close  to  the  point  of  solute  introduction, 
the  Impact  of  boundary  conditions  on  model  predictions  is  small. 

As  was  noted  in  Chapter  3,  the  assumption  that  no  mass  is  initially 
associated  with  the  immobile  region  has  a  large  effect  on  the  mobile 
region  plume  mass  and  velocity  behavior.  A  related  assumption  is  that 
water  which  is  obtained  from  the  sampling  wells  is  mobile  water.  The 
former  assumption  implies,  and  the  latter  states,  that  water  which  is 
injected  into  the  aquifer  by  the  injection  system,  or  withdrawn  from  the 


aquifer  by  the  sampling  system,  is  mobile  water.  Conceptually,  these 
assumptions  seem  reasonable.  If  it  is  assumed  that  regions  of  immobile 
water  do  indeed  exist  in  the  Borden  aquifer  (within  dead-end  pores  or 
zones  of  low  permeability,  for  example),  it  follows  that  during  injec¬ 
tion  or  extraction,  mobile  water  will  be  preferentially  displaced.  In 
light  of  this  preferential  displacement,  the  initial  condition  assump¬ 
tion  and  the  assumption  regarding  the  sampling  water  appear  reasonable. 
Therefore,  in  the  following  discussion,  the  moments  calculated  from  the 
spatial  data  will  be  assumed  to  be  mobile  region  plume  moments. 

Figures  5. A  through  5.7  compare  the  mass  in  solution  and  first 
spatial  moments  of  the  four  mobile  region  solute  plumes  (Freyberg,  1986; 
Roberts  et  al.,  1986)  with  the  moments  simulated  using  the  first-order 
rate  model.  The  Table  5.3  parameter  values  were  used  in  the  model. 
Also  shown  in  Figures  5. A  through  5.7  are  the  moments  simulated  using 
the  local  equilibrium  model.  The  local  equilibrium  model  simulations 
use  the  values  of  vQ,  ,  p,  6,  and  from  Table  5.3,  and  implicitly 
assume  $  -  f-  1  (i.e.,  all  water  is  mobile  water).  As  these  figures 
show,  the  local  equilibrium  model  and  the  physical  nonequilibrium  model 
produce  virtually  identical  results  for  the  parameters  used.  These 
results  indicate  that,  based  on  the  assumed  diffusion  coefficients  and 
grain  sizes,  for  the  time  scale  considered,  local  equilibrium  is  a  valid 
assumption.  In  particular,  notice  that  varying  the  first-order  rate 
constant  over  two  orders  of  magnitude  has  no  effect  on  the  results.  The 
reason  for  this  insensitivity  to  the  rate  constant  is  that  even  for  the 
minimum  value  of  the  rate  constant,  local  equilibrium  is  valid.  There¬ 
fore,  increasing  the  rate  constant  has  no  effect  on  the  simulations. 
Note,  though,  that  the  qualitative  fit  of  the  models  to  the  data  is 
poor.  The  mass  loss  and  the  deceleration  of  the  solute  plumes  which  was 
observed  experimentally  is  not  simulated  by  the  models.  One  way  to 
simulate  the  observed  spatial  moment  behavior  is  to  adjust  the  diffusion 
rate  constant.  Either  a  tenfold  increase  in  the  maximum  median  grain 
size  (from  b  ■  0.35  mm  to  b  ■  3.5  mm,  while  holding  constant)  or 
a  one-hundredfold  decrease  in  the  effective  intraparticle  diffusion 
coefficient  (while  holding  b  ■  0.35  mm  constant)  will  lead  to  a  one¬ 
hundredfold  decrease  in  the  minimum  rate  constant.  Figures  5.8  through 
5.11  compare  the  simulations  of  the  first-order  rate  model  using 
this  decreased  rate  constant  with  the  experimental  observations.  As  the 
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Figure  5.4.  Comparison  of  the  bromide  mobile  region  plume:  a)  Mass  in 
solution,  end  b)  First  spatial  moment,  estimated  from  syn¬ 
optic  data  and  simulated  using  an  equilibrium  and  a  first- 
order  rate  model  equivalent  to  a  spherical  diffusion  model. 
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Figure  5.5.  Comparison  of  the  bromoform  mobile  region  plume:  a)  Mass 
in  solution,  and  b)  First  spatial  moment,  estimated  from 
synoptic  data  and  simulated  using  an  equilibrium  and  a 
first-order  rate  model  equivalent  to  a  spherical  diffusion 
model. 
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Figure  5.6. 


Comparison  of  the  carbon  tetrachloride  mobile  region  plume: 
a)  Mass  in  solution,  and  b)  First  spatial  moment,  estimated 
from  synoptic  data  and  simulated  using  an  equilibrium  and  a 
firat-order  rate  model  equivalent  to  a  spherical  diffusion 
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Figure  5.7. 


Comparison  of  the  tetrachloroethylene  mobile  region  plume: 
a)  Mass  in  solution  and  b)  First  spatial  moment,  estimated 
from  synoptic  data  and  simulated  using  an  equilibrium  and  a 
first-order  rate  model  equivalent  to  a  spherical  diffusion 
model. 
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□  Estlmoted  from  data 
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Figure  5.9.  Comparison  of  the  bromoform  mobile  region  plume:  a)  Mass 
in  solution,  and  b)  First  spatial  moment,  estimated  from 
synoptic  data  and  simulated  using  a  first-order  rate  model 
equivalent  to  a  spherical  diffusion  model,  with  a  very  low 
rate  constant  (o'  ■  0.00594  d"*). 
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Figure  5.10.  Comparison  of  the  carbon  tetrachloride  mobile  region  plume: 

a)  Mass  in  solution,  and  b)  First  spatial  moment,  esti¬ 
mated  from  synoptic  data  and  simulated  using  a  first-order 
rate  model  equivalent  to  a  spherical  diffusion  model,  with 
a  very  low  rate  constant  (o'  ■  0.00608  d"*). 
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Figure  5.11.  Comparison  of  the  tetrachloroethylene  mobile  region  plume: 

a)  Mass  In  solution,  and  b)  First  spatial  moment,  esti¬ 
mated  from  synoptic  data  and  simulated  using  a  first-order 
rate  model  equivalent  to  a  spherical  diffusion  model,  with 
a  very  low  rate  constant  (o'  •  0.00566  d"*). 
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figures  show,  the  salient  characteristics  of  the  data  are  simulated  by 
the  model.  The  mass  in  mobile  solution  moment  simulations  demonstrate 
good  fits  to  the  data,  except  in  the  case  of  bromoform.  The  model, 
which  predicts  no  mass  loss  of  bromide  over  the  course  of  the  experi¬ 
ment,  is  consistent  with  the  data.  The  mass  loss  of  carbon  tetrachlo¬ 
ride  and  tetrachloroethylene  is  simulated  by  the  model.  The  observed 
mass  loss  of  bromoform,  however,  is  significantly  greater  than  model 
predictions.  This  discrepancy  may  be  due  to  the  biotransformation  of 
bromoform,  a  possibility  which  was  discussed  by  Roberts  et  al.  (1986). 
The  observed  deceleration  of  the  organic  plumes,  and  the  constant  veloc¬ 
ity  of  the  bromide  plume  are  qualitatively  simulated. 

The  only  parameter  which  was  adjusted  to  fit  the  model  to  the  data 
in  Figures  5.8  through  5.11  was  the  rate  constant  describing  solute 
transfer  into  the  regions  of  immobile  water.  The  fitted  value  of  this 
rate  constant  is  lOOx  less  than  was  calculated  using  literature  correla¬ 
tions  and  the  maximum  median  grain  size.  However,  this  hundredfold  dis¬ 
crepancy  is  not  entirely  unrealistic.  Hutzler  et  al.  (1986),  using  a 
similar  spherical  diffusion  model  to  simulate  soil  column  breakthrough 
responses,  found  that  grain  sizes  on  the  order  of  7-20x  greater  than  the 
mean  particle  grain  size  were  required  to  get  a  good  fit  of  the  model 
simulations  to  the  experimental  data.  This  Increase  in  grain  size  is 
equivalent  to  a  49-  to  400-fold  decrease  in  the  rate  constants  calcu¬ 
lated  using  literature  correlations.  An  alternate  explanation  is  that 
strongly  hindered  diffusion  is  occurring.  Such  hindered  diffusion  has 
been  reported  in  the  literature.  Roberts  (1966)  estimated  that  hexane 
diffusivlty  in  zeolite  crystals  was  seven  orders  of  magnitude  less  than 
bulk  diffusivlty.  Satterfield  et  al.  (1973)  found  that,  for  solute 
molecular  diameters  on  the  same  order  as  sorbent  pore  diameters,  liquid- 
phase  diffusivitles  of  hydrocarbons  were  less  than  bulk  diffusivlty  val¬ 
ues  by  several  orders  of  magnitude.  Prasher  and  Ma  (1977)  cited  other 
studies  where  the  reported  diffusivlty  values  for  organic  solutes  in  zeo¬ 
lites  were  several  orders  of  magnitude  less  than  bulk  diffusivlty  values. 
In  addition.  Ball  and  Roberts'  (1985)  finding  that  the  attainment  of 
equilibrium  in  batch  sorption  studies  measuring  tetrachloroethylene 
uptake  by  Borden  aquifer  material  required  at  least  several  months  also 
supports  the  hypothesis  of  strongly  hindered  diffusion  within  Borden 
aquifer  material. 
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Layered  Diffusion  Model 

Based  on  Sudicky's  (1986)  finding  of  low  permeability  horizontal 
lenses  within  the  Borden  aquifer,  it  is  also  reasonable  to  suppose  that 
solute  diffusion  into  these  layers  may  influence  transport.  Table  5.4 
lists  parameter  values  for  use  in  a  layered  diffusion  model.  The  ratio 
of  mobile  to  total  water  content  ($)  was  obtained  from  a  semi-independent 
measurement.  Freyberg  (1986)  showed  that  the  velocity  of  the  conserva¬ 
tive  tracers  averaged  0.091  m/d.  The  value  of  vQ  calculated  based  on 
hydraulic  conductivity  and  hydraulic  gradient  estimates,  and  assuming  a 
porosity  of  0.33,  was  0.079  m/d  (Mackay  et  al.,  1986).  If  it  is  assumed 
that  the  difference  between  the  two  velocity  estimates  is  caused  by 
zones  of  Immobile  water  which  do  not  contribute  to  transport,  the  value 
of  $  can  be  calculated  as  follows: 

_  0.079  m/d  _ 

*  0T091  m7d  0,86 


TABLE  5.4 

PARAMETER  VALUES  FOR  USE  IN  A  LAYERED  DIFFUSION  MODEL 


Bulk  soil  density  (p) 
Total  porosity  (0) 

Ratio  of  mobile  to  total 

water  content 

(♦> 

1.81  g/cm3 

0.33 

0.86 

Pore  water  velocity  (vQ) 

Fraction  of  "mobile"  sorption  sites  (f) 

0.079  m/d 

0.86 

Layer  half-width  (b) 

Bromide 

Bromo- 

form 

0.01-0.05  m 

Carbon  Tetra- 

Tetra-  chloro- 

chloride  ethylene 

Distribution  coeff. 

(Kd),  cm3/g 

0.0 

0.17 

0.17  0.48 

Mass  of  injected 
solute  (Mj) ,  g 

3,870 

0.38 

0.37  0.36 

Liquid  diffusion 
coeff.  (DQ),  m2/d 

1.74xl0”4 

7.17xlO”5 

7.34xl0”5  6.83xl0"5 

Effective  diffusion 
coeff.  (Dg) ,  nr/d 

2. 9-8. 7 

1.2-3 .6 

1. 2-3.7  1.1-3. 4 

xlO"5 

xlO”5 

xlO”5  xlO"5 

Equivalent  first-order 
rate  constant  (o'). 

1.57xl0-3 

6.51xl0”4 

6.63xl0"4  6.18xl0"4 

d"1 

-  1. 18x10” 1 

-  4.86xl0“2 

-  4.98xl0”2  -  4.65xl0"2 
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Of  course,  It  oust  be  realized  that  the  value  of  vQ  represents  a  rough 
estimate,  so  that  the  value  of  $  calculated  above  should  be  viewed  as  a 
gross  approximation.  However,  for  this  exercise,  which  is  meant  to  com¬ 
pare  qualitative  characteristics  of  model  simulations  with  the  data, 
such  an  approximation  is  adequate. 

In  the  case  of  the  layered  model,  diffusion  is  not  assumed  to  be 
intraparticle  as  it  was  for  the  spherical  model,  but  rather  interparti- 
cle.  For  interparticle  diffusion,  values  of  x  2-6  have  been  reported 
(Cussler,  1984).  Thus,  a  range  of  may  be  calculated  by  dividing  DQ 
for  the  various  compounds  by  a  factor  of  2  to  6.  The  range  of  immobile 
layer  half-widths  (b)  was  estimated  based  on  the  approximate  dimensions 
of  the  low  permeability  layers  noted  in  Freyberg  (1986) .  The  fraction 
of  sorption  sites  adjacent  to  zones  of  mobile  water  (f)  was  approximated 
assuming  an  equal  distribution  of  sorption  sites  within  the  mobile  and 
immobile  region  volumes.  Other  parameter  values  are  the  same  as  those 
used  with  the  spherical  diffusion  model. 

As  with  the  spherical  diffusion  model,  it  is  now  possible  to  calcu¬ 
late  first-order  rate  constants  (o')  equivalent  to  the  diffusion  model 
rate  constants  which  can  be  obtained  from  the  Table  5.4  parameter 
values.  The  range  in  the  values  of  o'  is  due  to  the  range  in  the  values 
of  Dg  and  b,  with  the  maximum  a’  calculated  using  the  maximum  value  of 
Dg  and  the  minimum  value  of  b,  and  the  minimum  a’  obtained  using  the 
minimum  and  the  maximum  b. 

Figures  5.12  through  5.15  compare  the  mass  in  solution  and  first 
spatial  moments  calculated  from  the  data  of  the  four  mobile  region  sol¬ 
ute  plumes  with  the  mass  in  solution  and  first  moments  simulated  using 
the  first-order  rate  model  equivalent  to  the  layered  diffusion  model. 
Also  shown  are  simulations  using  the  local  equilibrium  model.  As  in 
Figures  5.4  through  5.7  depicting  the  spherical  diffusion  model,  the 
differences  between  the  nonequilibrium  and  the  equilibrium  simulations 
do  not  appear  significant,  though  at  least  for  the  low  rate  constant 
simulations,  the  salient  aspects  of  the  solute  behavior  (mass  loss  and 
deceleration)  are  seen.  The  magnitude  of  the  simulated  behavior,  how¬ 
ever,  is  considerably  less  than  that  exhibited  by  the  actual  plumes. 
The  insignificance  of  the  effect  of  the  immobile  region  is  largely  due 
to  the  assumption  of  the  equal  distribution  of  sorption  sites  between  the 
mobile  and  immobile  water  regions.  Due  to  this  assumption,  the  effect  of 
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Figure  5.12. 


Comparison  of  the  bromide  mobile  region  plume:  a)  Mass  in 
solution,  and  b)  First  spatial  moment,  estimated  from  syn¬ 
optic  data  and  simulated  using  an  equilibrium  and  a  first- 
order  rate  model  equivalent  to  a  layered  diffusion  model. 
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Figure  5.13.  Comparison  of  the  bromoform  mobile  region  plume:  a)  Mass 
in  solution,  and  b)  First  spatial  moment,  estimated  from 
synoptic  data  and  simulated  using  an  equilibrium  and  a 
first-order  rate  model  equivalent  to  a  layered  diffusion 
model. 
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Figure  5.14.  Comparison  of  the  carbon  tetrachloride  mobile  region  plume: 

a)  Mass  in  solution  and  b)  First  spatial  moment,  estimated 
from  synoptic  data  and  simulated  using  an  equilibrium  and 
a  first-order  rate  model  equivalent  to  a  layered  diffusion 
model. 
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Comparison  of  the  tetrachloroethylene  mobile  region  plume: 
a)  Mass  in  solution  and  b)  First  spatial  moment,  estimated 
from  synoptic  data  and  simulated  using  an  equilibrium  and 
a  first-order  rate  model  equivalent  to  a  layered  diffusion 
model. 


the  immobile  region  on  solute  transport  -is  limited,  since  relatively 
little  solute  is  associated  with  the  immobile  regions.  Figures  5.16 
through  5.19  demonstrate  the  impact  of  decreasing  the  value  of  f,  thus 
implying  that  more  sorption  sites  are  within  regions  of  low  permeabil¬ 
ity.  The  assumption  that  low  permeability  regions,  perhaps  made  up  of 
silty  material,  have  a  higher  sorption  capacity  is  not  unreasonable, 
though  at  the  Borden  site,  there  is  no  evidence  correlating  low  hydrau¬ 
lic  conductivity  with  high  distribution  coefficients  (Durant,  1986). 

As  Figures  5.16  through  5.19  show  (using  a  value  of  f  *  0.4  and  a 
value  of  the  rate  constant  based  on  average  values  of  Dg  and  b),  the 
main  characteristics  of  the  data  can  be  simulated,  again  with  the  excep¬ 
tion  of  the  decrease  in  bromoform  mass  (Figure  5.17),  which  is  attrib¬ 
uted  to  transformation. 

If  diffusion  into  low  permeability  lenses  is  occurring,  and  if 
advective  transport  through  the  lenses  is  faster  than  diffusive  trans¬ 
port  into  and  out  of  the  lenses,  then  advective  transport  would  dominate 
and  the  effect  of  diffusion  into  the  lenses  would  be  minor. 

The  concept  of  comparing  advective  and  diffusive  rates.  Introduced 
in  Chapter  3,  may  be  used  to  determine  the  relative  Importance  of  advec¬ 
tive  and  diffusive  transport  into  low  permeability  lenses.  Assume 
lenses  that  are  0.06  m  thick  (b  -  0.03  m),  with  hydraulic  conductivity 
two  orders  of  magnitude  less  than  the  hydraulic  conductivity  of  the 
"mobile"  regions.  Therefore: 


mm  9-i  *  l0"A 


and  the  advective  rate  constant  for  tetrachloroethylene  transport  in  the 
lenses  would  be: 


vlens  _  (9.1x10  *  m/d)  _  ,  w  ,rt-3  j-1 

2b  R,  2(0.03  m)(12. 5)  1,21  10  d 

lm 


The  diffusive  rate  constant  would  be: 


•»n»  _s  ? 

.  e  _  3(1.71x10  3  n>  /d) 

Rlmb2  (12.5)(0.03  m)2 
using  x  ■  4  to  calculate  from  D0. 
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Figure  5.16.  Comparison  of  the  bromide  mobile  region  plume:  a)  Mass  In 
solution,  and  b)  First  spatial  moment,  estimated  from  syn¬ 
optic  data  and  simulated  using  a  first-order  rate  model 
equivalent  to  a  layered  diffusion  model,  with  high  sorp¬ 
tion  capacity  within  the  immobile  regions  (f  •  0.4). 
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Figure  5.17.  Comparison  of  the  bromoform  mobile  region  plume:  a)  Maas 
in  solution,  and  b)  First  spatial  moment,  estimated  from 
synoptic  data  and  simulated  using  s  first-order  rate  model 
equivalent  to  a  layered  diffusion  model,  vith  high  sorp¬ 
tion  capacity  within  the  immobile  regions  (f  *  0.4). 


95 


0.35 


□  Estimated  from  data 


Increased  immobile  sorption  site  simulation 


</>  0.2 


TIME  (days) 


40  -J  □  Estimated  from  data 


Increased  immobile  sorption  site  simulation 


TIME  (days) 


Figure  5.18.  Comparison  of  the  carbon  tetrachloride  mobile  region  plume: 

a)  Mass  in  solution,  and  b)  First  spatial  moment,  estimated 
from  synoptic  data  and  simulated  using  a  first-order  rate 
model  equivalent  to  a  layered  diffusion  model,  with  high 
sorption  capacity  within  the  immobile  regions  (f  “  0.4). 
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Figure  5.19.  Comparison  of  the  tetrachloroethylene  mobile  region  plume: 

a)  Mass  in  solution,  and  b)  First  spatial  moment,  estimated 
from  synoptic  data  and  simulated  using  a  first-order  rate 
model  equivalent  to  a  layered  diffusion  model,  with  high 
sorption  capacity  within  the  Immobile  regions  (f  *  0.4). 


Comparing  the  tvo  rate  constants,  it  is  seen  that  for  the  assumed 
parameter  values,  the  diffusive  rate  constant  is  nearly  four  times 
greater  than  the  advective  rate  constant.  This  indicates  that  for  rea¬ 
sonable  assumptions  regarding  lens  characteristics  (b  >_  0.03  m;  lens 
hydraulic  conductivity  <_ 0.01  mobile  region  hydraulic  conductivity),  the 
effect  of  diffusive  transport  would  dominate  the  effect  of  advective 
transport  within  the  lenses.  As  the  ratio  of  advective  to  diffusive 
rate  constants  is  independent  of  R^n,  the  dominance  of  diffusive  trans¬ 
port  applies  for  all  the  compounds,  since  of  all  the  solutes,  tetra- 
chloroethylene  has  the  lowest  estimated  value  of  D^. 


Temporal  Data  Simulation 

Another  way  of  demonstrating  the  complementarity  of  spatial  and 
temporal  transport  behavior  is  to  use  the  model  parameters  obtained  in 
the  preceding  analysis  of  the  zeroth  and  first  spatial  moments  to  simu¬ 
late  temporal  responses.  The  following  discussion  will  use  parameters 
obtained  assuming  that  layered  diffusion  dominates  transport  (Table  5.4 
with  f  ■  0.40),  although  either  the  layered  or  spherical  diffusion  model 
assumptions  could  have  been  used  for  this  exercise. 

Values  for  longitudinal  and  transverse  dispersion  coefficients  and 

Initial  plume  dimensions  must  be  determined  for  use  as  model  input 

parameters.  These  parameters  may  be  estimated  from  the  second  spatial 

2  2 

moment  data.  Figures  5.20  through  5.27  show  o^  and  o^y,  the  principal 
values  of  the  spatial  covariance  tensor,  calculated  for  the  four  solutes 
at  each  synoptic  sampling  time.  Figures  5.20a  through  5.27a  compare 
simulations  of  the  local  equilibrium  model  with  principal  values  of  the 
covariance  tensor.  The  values  of  the  longitudinal  and  transverse  dis¬ 
persion  coefficients  (D^  and  Dy)  and  the  Initial  plume  length  and  width 
(2L  and  2M)  were  chosen  to  provide  simultaneously  a  visual  best  fit  to 
the  eight  data  sets.  The  other  parameter  values  used  in  the  equilibrium 
sndel  (vQ,  K^,  p  and  9)  were  obtained  in  the  preceding  section.  Simi¬ 
larly,  Figures  5.20b  through  5.27b  compare  the  second  spatial  moment  data 
with  simulations  of  the  first-order  rate  model.  Recall  that  the  first- 
order  rate  constant  used  in  the  model  was  calculated  to  be  equivalent  to 
the  layered  diffusion  model  rate  parameter.  As  with  the  local  equilib¬ 
rium  model,  the  longitudinal  and  transverse  dispersion  coefficients  in 
the  mobile  region  (Dw  and  D0y)  and  the  initial  plume  dimensions  (2L  and 
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Figure  5.20.  Comparison  of  the  bromide  mobile  region  plume  principal 
component  of  the  spatial  covariance  tensor  in  the  longi¬ 
tudinal  direction  (c^)  estimated  from  the  synoptic  data 
and  simulated  by  fitting  the  a)  Equilibrium  model,  and 
b)  First-order  rate  model  to  the  estimates. 
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Figure  5.21.  Comparison  of  the  bromoform  mobile  region  plume  principal 
component  of  the  spatial  covariance  tensor  in  the  longi¬ 
tudinal  direction  (oxx)  estimated  from  the  synoptic  data 
and  simulated  by  fitting  the  a)  Equilibrium  model,  and 
b)  First-order  rate  model  to  the  estimates. 
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Figure  5.22.  Comparison  of  the  carbon  tetrachloride  mobile  region  plume 
principal  component  of  the  spatial  covariance  tensor  in  the 
longitudinal  direction  (o^)  estimated  from  the  synoptic 
data  and  simulated  by  fitting  the  a)  Equilibrium  model, 
and  b)  First-order  rate  model  to  the  estimates. 
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Figure  5.23. 


Comparison  of  the  tetrachloroethylene  mobile  region  plume 
principal  component  of  the  spatial  covariance  tensor  in  the 
longitudinal  direction  (o^)  estimated  from  the  synoptic 
data  and  simulated  by  fitting  the  a)  Equilibrium  model, 
and  b)  First-order  rate  model  to  the  estimates. 
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Figure  5.24.  Comparison  of  the  bromide  mobile  region  plume  principal 
component  of  the  spatial  covariance  tensor  in  the  trans¬ 
verse  direction  (°yy)  estimated  from  the  synoptic  data 
and  simulated  by  f feting  the  a)  Equilibrium  model,  and 
b)  First-order  rate  model  to  the  estimates. 
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Figure  5.25.  Comparison  of  the  bromoform  mobile  region  plume  principal 
component  of  the  spatial  covariance  tensor  in  the  trans¬ 
verse  direction  (o_y)  estimated  from  the  synoptic  data 
and  simulated  by  fitting  the  a)  Equilibrium  model,  and 
b)  First-order  rate  model  to  the  estimates. 
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Figure  5.26.  Comparison  of  the  carbon  tetrachloride  mobile  region  plume 
principal  component  of  the  spatial  covariance  tensor  in  the 
transverse  direction  (o^  )  estimated  from  the  synoptic 
data  and  simulated  by  fitting  the  a)  Equilibrium  model, 
and  b)  First-order  rate  model  to  the  estimates. 
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Figure  5.27. 


Comparison  of  Che  tetrachloroethylene  mobile  region  plume 
principal  component  of  Che  spatial  covariance  tensor  in  the 
transverse  direction  (°yy)  estimated  from  the  synoptic 
data  and  simulated  by  fitting  the  a)  Equilibrium  model, 
and  b)  First-order  rate  model  to  the  estimates. 


2M)  were  adjusted  to  obtain  simultaneously  a  visual  best  fit  to  the 
eight  data  sets  being  considered.  Other  parameters  used  in  the  model 
were  obtained  in  the  previous  section  (Table  5. A).  A  value  of  f  *  0.4 
was  used  to  test  whether  parameters  which  have  been  shown  to  simulate 
the  spatial  behavior  can  also  simulate  the  temporal  behavior.  Table  5.5 
lists  the  values  for  the  dispersion  coefficients  and  the  initial  plume 
dimensions  which  will  be  used  in  the  following  analysis  of  the  temporal 
behavior. 


TABLE  5.5 


VALUES  FOR  DISPERSION  COEFFICIENTS  AND  INITIAL  PLUME 
DIMENSIONS  OBTAINED  FROM  SPATIAL  SECOND  MOMENT  DATA 


Local  Equilibrium  Model 

D;  -  0.04  m2/d 

D'  -  0.008  m2/d 

L  -  2.5  m 

M  -  2.0  m 

First-Order  Rate  Model* 

Dm„  -  0.02  m2/d 
mx  a 

Dmy  ■  0.008  nr/d 

L  -  2.5  m 

M  ■  2.0  m 

*Equivalent  to  layered  diffusion  model. 

Using  the  parameter  values  in  Tables  5.4  and  5.5,  it  is  possible  to 
simulate  both  the  temporal  moments  and  the  measured  temporal  responses 
at  particular  sampling  wells.  Figures  5.28  through  5.31  compare  zeroth 
and  first  temporal  moments  estimated  from  the  data  with  predictions  of 
equilibrium  and  nonequilibrium  models.  Moment  estimates  were  obtained 
by  numerically  evaluating  the  breakthrough  response  data  at  the  various 
sampling  wells,  depicted  in  Appendix  E,  using  the  formulae: 

n 


o,t 


3-1 


CJ(t3  "  tj~l) 


(5-3) 


I 


Vj(tj 


) 


where  Cj  *  the  normalized  concentration  measured  at  time  tj. 
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Figure  5.28.  Comparison  of  bromide  temporal  moments  estimated  from  the 
data  and  predicted  by  equilibrium  and  nonequilibrium 
models:  a)  Zeroth,  and  b)  First  moments. 
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Figure  5.29.  Comparison  of  bromoform  temporal  moments  estimated  from 
the  data  and  predicted  by  equilibrium  and  nonequilibrium 
models:  a)  Zeroth,  and  b)  First  moments. 
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Figure  5.30.  Comparison  of  carbon  tetrachloride  temporal  moments  esti 
mated  from  the  data  and  predicted  by  equilibrium  and  non 
equilibrium  models:  a)  Zeroth,  and  b)  First  moments. 
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Figure  5.31.  Comparison  of  cetrachloroethylene  temporal  moments  esti' 
mated  from  the  data  and  predicted  by  equilibrium  and  non' 
equilibrium  models:  a)  Zeroth,  and  b)  First  moments. 
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Model  predictions  were  aade  using  the  equilibrium  and  nonequilib¬ 
rium  models  to  simulate  breakthrough  responses*  The  moments  of  these 
responses  were  then  evaluated  using  the  same  numerical  technique  that 
had  been  applied  to  the  data.  In  Chapter  3,  Table  3.3,  it  was  shown 
that  both  equilibrium  and  nonequilibrium  models  have  the  same  expression 
for  the  zeroth  temporal  moment  in  terms  of  the  model's  parameters.  Thus, 
differences  in  the  zeroth  moment  simulations  of  the  two  model  types  are 
due  to  parameter  differences,  and  not  to  model  structure.  Although 
certain  trends  are  evident,  the  zeroth  moment  data  and  simulations  show 
a  fair  amount  of  scatter.  Both  models  simulate  the  low  values  of  the 
zeroth  moment  found  at  the  sampling  well  located  at  x  *  4.79  m.  This 
low  value  is  observed  because  the  sampling  well  is  offset  a  considerable 
distance  from  the  center  line  of  the  plume.  It  is  also  apparent  that 
the  zeroth  moment  data  values  for  bromoform  at  the  far-fleld  wells  are 
considerably  lower  than  the  simulated  values.  This  discrepancy  may  be 
due  to  the  possible  biotransformation  of  bromoform,  which  was  discussed 
earlier  in  connection  with  the  zeroth  spatial  moment  results. 

From  Table  3.3  in  Chapter  3,  it  was  also  found  that  the  first  tem¬ 
poral  moment  calculated  using  the  nonequilibrium  models  is  equal  to  the 
first  temporal  moment  of  the  equilibrium  model  multiplied  by  a  constant. 
Again,  therefore,  any  differences  in  the  first  temporal  moment  simula¬ 
tions  are  due  to  parameter  differences,  rather  than  to  differences  in 
model  structure.  Figures  5.28b  through  5.31b  show  the  first  temporal 
moment  as  a  function  of  distance  in  the  x  direction.  Reasonably  good 
agreement  is  found  between  simulated  moments  and  moments  calculated  from 
the  data. 

Judging  from  the  zeroth  and  first  temporal  moment  simulations, 
there  is  little  evidence  for  saying  one  model  or  the  other  is  better  at 
describing  the  data,  largely  because  both  type  models  can  be  considered 
"pseudo-equilibrium  models"  with  respect  to  temporal  moments.  That  is, 
since  the  zeroth  and  first  temporal  moments  calculated  using  the  non- 
equlllbrlum  models  are  independent  of  the  rate  parameter,  the  zeroth  and 
first  moments  can  be  reproduced  by  an  equivalent  equilibrium  model. 
Figures  5.28  through  5.31  do  indicate,  ho”ever,  that  the  parameter  val¬ 
ues  obtained  from  fitting  the  spatial  moment  data  can  be  used  to  predict 
the  temporal  moment  data  fairly  well,  for  both  equilibrium  and  nonequi¬ 
librium  models. 
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Another  way  of  examining  the  temporal  data  is  to  inspect  the  break¬ 
through  responses  themselves,  at  various  wells.  Using  one-dimensional 
models,  Goltz  and  Roberts  (1986)  showed  that  the  sharp  breakthrough  and 
tailing  exhibited  by  the  near-field  well  breakthrough  response  data  was 
better  simulated  using  nonequilibrium  than  equilibrium  models.  Figures 
5.32  through  5.39  compare  simulations  of  three-dimensional  equilibrium 
and  first-order  rate  models  with  the  data.  For  each  solute,  two  wells 
were  chosen  for  comparison,  one  near-field  well  and  one  far-field  well. 
Figure  5.37,  which  compares  simulations  of  the  equilibrium  and  first- 
order  rate  models  with  the  near-field  well  response  data  for  bromoform, 
also  includes  a  layered  diffusion  model  simulation.  The  similarity  of 
the  first-order  rate  and  layered  diffusion  model  simulations  provides 
further  evidence  justifying  the  use  of  the  first-order  rate  model  to 
approximate  the  more  complex  diffusion  models. 

As  the  temporal  moments  have  already  been  discussed,  the  following 
analysis  will  concentrate  on  the  general  characteristics  and  form  of  the 
breakthrough  responses.  The  breakthrough  responses  simulated  for  the 
four  solutes  at  the  far-field  wells  are  fairly  symmetric,  with  the  non- 
equilibrium  model  response  exhibiting  a  bit  more  asymmetry  than  that  of 
the  equilibrium  model.  It  is  interesting  to  note  that  for  the  bromide 
simulated  responses,  the  equilibrium  model  curve  exhibits  more  spreading 
and  a  lower  peak  value  than  that  of  the  nonequilibrium  model,  whereas 
for  the  organic  simulations,  the  situation  is  reversed.  The  reason  for 
this  behavior  may  be  found  by  applying  Eq.  3-66  to  the  nonequilibrium 
model  parameter  values  for  the  different  solutes.  Applying  Eq.  3-66, 
which  calculates  an  equivalent  dispersion  coefficient  (Deff),  reveals 
Deff  for  the  bromide  at  the  sampling  well  to  be  less  than  Deff  for  the 
organics.  The  value  of  Deff  which  is  calculated  using  Eq.  3-66  is  an 
equilibrium  model  equivalent  parameter  which  combines  the  effect  of 
dispersion  in  the  mobile  region  with  the  effect  of  so-called  "holdup 
dispersion”  (Koch  and  Brady,  1986).  (Holdup  dispersion  describes  a 
spreading  mechanism  caused  by  solute  diffusion  into  immobile  regions.) 
Owing  to  the  sorption  of  the  organic  solutes  within  the  immobile  region, 
the  impact  of  holdup  dispersion  on  the  organic  compounds'  behavior  is 
more  pronounced  than  the  impact  on  the  conservative  tracer's  behavior, 
so  the  organic  simulations  exhibit  relatively  more  spreading  at  the  far- 
field  sampling  points. 
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Figure  5.32. 


Breakthrough  response  data  and  model  predictions  for 
bromide  at  a  far-field  well  (x  *  21.0  m,  y  ■  9.0  m, 
z  *  -A. 17  m). 
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Figure  5.33. 


Breakthrough  response  data  and  model  predictions  for 
bromoform  at  a  far-field  well  (x  ”  21.0  m,  y  *  9.0  m, 
z  ■  -A. 17  m). 
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NORMALIZED  CONCENTRATION 


Examination  of  the  near-field  responses  seems  to  confirm  Goltz  and 
Roberts'  (1986)  conclusion  that  the  characteristically  sharp  break¬ 
through  and  tailing  of  the  breakthrough  response  data  are  better  simu¬ 
lated  with  the  nonequilibrium  models.  On  the  other  hand,  in  the  far 
field  neither  model  type  offers  a  significant  advantage  over  the  other 
for  simulating  the  data.  This  is  to  be  expected,  since  at  long  dis¬ 
tances,  the  nonequilibrium  models  can  be  approximated  by  equivalent 
equilibrium  models. 

Overall,  considering  the  fact  that  input  parameters  were  not  fitted 
to  the  temporal  data,  both  type  models  quite  adequately  simulate  the 
observed  data.  However,  especially  in  the  case  of  tetrachloroethylene , 
the  solute  expected  to  be  most  impacted  by  nonequilibrium  effects  due  to 
its  high  distribution  coefficient,  the  salient  characteristics  of  the 
near-field  data  seem  to  be  better  captured  by  the  nonequilibrium  model 
simulations.  Also  Figure  5.37  shows  little  difference  between  first- 
order  rate  and  diffusion  model  breakthrough  simulations,  certainly  not 


enough  difference  to  choose  one  model  over  the  other  based  solely  on  the 


goodness  of  fit. 


ALTERNATE  HYPOTHESES 


Although  the  hypothesis  of  physical  nonequilibrium  appears  to 
explain  tne  observed  temporal  and  spatial  behavior  of  the  sorbing  and 
nonsorbing  solutes,  there  are  other  hypotheses  which  also  may  adequately 
explain  the  observations.  This  section  examines  some  of  these  alterna¬ 
tive  hypotheses,  in  light  of  what  is  known  about  the  Borden  situation. 


Nonlinear  Sorption 


If  solute  transport  is  assumed  to  be  governed  by  the  advective /dis¬ 
persive  equation,  with  sorption  described  by  a  Freundlich  Isotherm  (van 
Genuchten  and  Cleary,  1982),  the  following  equations,  written  for  one- 
dimensional  transport,  apply: 


3C  .  n,  32C  3C  p  3S 

3t  ux  TT  o  3x  “  9  3t 


(5-4a) 


S  -  KCn 


(5-4b) 


For  n  <  1,  sharp  breakthrough  and  tailing  are  characteristic  of  the 
breakthrough  responses,  making  them  similar  to  those  seen  at  Borden  (van 
Genuchten  and  Cleary,  1982).  These  breakthrough  response  characteristics 
are  due  to  the  fact  that  at  high  aqueous  concentrations,  the  sorption 
partition  coefficient  is  smaller  than  at  low  concentrations,  so  that  the 
main  body  of  the  solute  pulse  overtakes  the  leading  edge.  On  the  other 
hand,  the  low  concentrations  in  the  trailing  edge  of  the  pulse  are  sorbed 
more  strongly  than  the  main  body  of  the  pulse,  creating  a  long  tail. 

Qualitatively,  at  least,  this  explanation  accounts  for  the  relative 
symmetry  of  the  nonsorbing  bromide  tracer  breakthrough  responses,  as 
well  as  the  sharp  breakthrough  and  tailing  of  the  organic  solutes' 
breakthrough  responses.  It  may  also  explain  the  loss  in  mass  and  the 
deceleration  of  the  organic  solutes;  as  the  solute  plume  disperses  over 
time,  the  aqueous  concentrations  decline,  and  the  overall  sorption 
increases.  To  examine  these  qualitative  statements  in  more  detail,  a 
numerical  analysis  was  performed. 

A  Crank -Nicolson  finite  difference  representation  of  Gq.  5-4  was 
developed.  Parameters  for  the  model  were  obtained  from  Curtis  et  al. 
(1986).  By  performing  batch  equilibrium  sorption  studies  using  the 
Borden  aquifer  material,  Curtis  et  al.  (1986)  found  that  at  the  low  con¬ 
centrations  seen  in  the  Borden  experiment,  the  sorption  of  the  organic 
compounds  could  be  modeled  by  a  linear  Isotherm.  However,  Curtis  et  al. 
(1986)  fit  their  isotherm  data  with  a  Freundlich  model  as  well.  They 
found  that  of  all  the  solutes,  tetrachloroethylene  had  a  Freundlich 
exponent  most  different  from  unity.  Using  the  following  parameter 
values  obtained  from  Curtis  et  all.'s  (1986)  sorption  experiments: 


K  -  0.93  (cm3/g)0*79 


n  -  0.79 


combined  with  the  following  hydrogeologic  parameter  values  (Tables  5.3- 
5.5): 

p  •  1.81  g/cm3 


6  -  0.33 


vQ  ■  0.079  m/d 


■  0.04  m3/d 


'•/'o  •' 


simulations  were  run  to  determine  the  zeroth  and  first  spatial  moments 
of  the  solute  plume  distribution  versus  sampling  time.  Figure  5.40 
compares  simulations  with  moments  estimated  from  the  tetrachloroethylene 
data.  The  decline  in  mass  and  deceleration  simulated  by  the  model  is 
much  less  than  exhibited  by  the  data.  The  reason  that  the  simulated 
loss  in  mass  and  deceleration  is  relatively  Insignificant  is  due  to  the 
assumption  of  equilibrium.  The  initial  solute  plume  is  assumed  to  be 
instantaneously  equilibrated,  with  a  fraction  of  the  solute  mass  imme¬ 
diately  sorbed.  Thus,  even  though  more  and  more  mass  goes  into  the 
sorbed  phase  with  time,  the  ratio  of  sorbed  to  aqueous  mass  does  not 
significantly  change. 


Hysteretic  Sorption 


Investigators  have  found  that  sorption  may  be  hysteretic  (Hornsby 
and  Davidson,  1973;  Swanson  and  Dutt,  1973;  van  Genuchten  et  al.  1974; 
Horzempa  and  DiToro,  1983;  DiToro,  1985).  Hysteretic  sorption/desorp¬ 
tion  isotherms,  combined  with  advective/dispersive  transport,  have  been 
shown  to  simulate  breakthrough  responses  with  long  tails  (van  Genuchten 
et  al.,  1974;  van  Genuchten  and  Cleary,  1982). 

To  Investigate  whether  sorption/desorption  hysteresis  is  an  impor¬ 
tant  mechanism  at  Borden,  Curtis  et  al.  (1986)  performed  a  short-term 
(6  day)  sorption/desorption  experiment  using  Borden  aquifer  material  and 
hexachloroethane.  No  significant  sorption/desorption  hysteresis  was 
observed . 

In  this  work,  a  long-term  (60  day)  sorption/desorption  study  was 
conducted  using  tetrachloroethylene  and  Borden  aquifer  material. 
Tetrachloroethylene  was  chosen,  as  it  is  the  most  strongly  sorbing  of 
the  solutes  under  consideration,  and  would  presumably  exhibit  the  most 
obvious  hysteresis.  A  long-term  study  was  conducted  based  on  the  possi¬ 
bility  of  slow  sorption  due  to  diffusion  rate  limitations.  Experimental 
details  are  presented  in  Appendix  F.  Figure  5.41  plots  the  sorption/de¬ 
sorption  isotherm  data. 

The  model  most  often  used  to  describe  hysteretic  sorption/desorp¬ 
tion  isotherms,  in  conjunction  with  one-dimensional  advective/dispersive 
solute  transport  (Hornsby  and  Davidson,  1973;  van  Genuchten  et  al., 
1974;  Wood  and  Davidson,  1975;  and  van  Genuchten  et  al.,  1977)  is: 
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Figure  5.40.  Comparison  of  the  tetrachloroethylene  plume  a)  mass  in 
solution,  and  b)  first  spatial  moment,  estimated  from  the 
data  and  simulated  by  an  equilibrium  transport  model 
assuming  nonlinear  sorption  and  using  experimentally 
obtained  sorption  parameter  values. 
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Figure  5.41.  Sorption/desorption  isotherm  results  for  uptake  of  tetra- 
chloroethylene  on  Borden  aquifer  material:  a)  Fit  using  a 
linear  sorptlon/hysteretic  desorption  isotherm.  b)  Fit 
using  a  nonlinear  sorption/hysteretic  desorption  isotherm. 


This  model  combines  the  advective/disperslve  solute  transport  equa¬ 
tion  with  a  sorption  source/slnk  term  desclbed  by  a  Freundlich  Isotherm. 
Depending  on  whether  solute  is  sorbing  or  desorbing,  different  sets  of 
Freundlich  isotherm  parameters  are  used. 

The  tetrachloroethylene  sorption/desorption  data  shown  In  Figure 
5.41  were  used  to  determine  Freundlich  isotherm  parameters  for  this 
model.  In  Figure  5.41a,  linear  sorption  is  assumed  (na(js  -  1.0)  and 
using  least-squares  regression,  best  fit  values  of  Ka(j8  and  ndeg  were 
determined.  In  Figure  5.41b,  both  sorption  and  desorption  are  described 
by  nonlinear  isotherms.  Least-squares  regression  was  used  to  determine 
best  fit  parameter  values  for  nads,  ndeg,  and  Kadg.  In  both  Figures 
5.41a  and  5.41b,  a  value  of  -  13.6  pg/1  was  assumed.  This  Cmax 
value  corresponds  to  the  maximum  aqueous  concentration  measured  in  the 
sorption  half  of  the  sorption/desorption  experiments. 


A  Crank-Nlcolson  finite  difference  approximation  of  Eq.  5-5  was 
developed.  Assuming  linear  sorption  and  hysteretlc  desorption  (Figure 
5.41a),  the  following  set  of  parameter  values  was  used  in  the  model: 


p  -  1.81  g/cm^ 


6  -  0.33 


Kads  "  °*69  cm  '  8 


Dx  »  0.04  m^/d 


nads  "  l-° 


vQ  ■  0.079  m/d 


ndes  "  0,66 


Assuming  nonlinear  sorption  and  hysteretlc  desorption  (Figure  5.41b), 
the  above  physical  and  hydrodynamic  parameters  were  used  with: 


Kads  -  1.18  <c«3/g)0-76 

nads  "  0,76 
ndes  "  °*59 


Figures  5.42  and  5.43  compare  model  simulations  of  the  zeroth  and 
first  spatial  moments  obtained  using  these  two  sets  of  parameter  values 
with  moments  estimated  from  the  data.  The  zeroth  and  first  moment  simu¬ 
lations  obtained  assuming  a  linear  sorption  isotherm  with  hysteretlc 
desorption  (Figure  5.42)  do  not  exhibit  the  observed  characteristics  of 
the  experimental  data.  The  mass  in  solution  exhibits  a  minimum  with 
respect  to  time,  and  the  first  moment  versus  time  plot  is  nearly  linear. 

The  zeroth  and  first  moment  simulations  obtained  using  the  non¬ 
linear  sorption  isotherm  with  hysteretlc  desorption  (Figure  5.43)  are 
qualitatively  comparable  to  the  observed  data,  though  quantitatively, 
the  simulated  loss  in  mass  and  deceleration  are  much  less  than  observed. 
The  qualitative  similarity,  however,  indicates  that  it  may  be  possible 


to  fit  such  a  nonlinear  sorption/hysteretlc  desorption  model  to  the 
data.  Figure  5.44  compares  model  simulations  with  moments  estimated 
from  the  data  for  the  following  Freundlich  isotherm  parameter  values 
which  were  selected  in  an  attempt  to  fit  the  data: 


Kads  "  0,50  (cm3/ g)^ 


-  0.70 


ndes  “  °*40 


As  Figure  5.44  shows,  the  nonlinear  sorption/hysteretlc  desorption  model 
can  adequately  simulate  the  main  characteristics  of  the  observed  spatial 


moment  data. 


Chemical  Nonequilibrium 


Chemical  nonequilibrium  models  that  describe  solute  transport  have 
been  proposed  (Lindstrom  and  Narasimhan,  1973;  Nkedi-Kizza  et  al.,  1984; 
Valocchl,  1985a).  One  such  model  combines  the  advective/dispersive 
transport  equation  with  a  first-order  rate  expression  to  describe  a 


chemical  reaction: 
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Figure  5.42.  Comparison  of  Che  tetrachloroethylene  plume  a)  mass  In 
solution,  and  b)  first  spatial  moment,  estimated  from  the 
data  and  simulated  by  an  equilibrium  transport  model  assum' 
ing  linear  sorption/hysteretic  desorption  and  using  exper-* 
lmentally  obtained  sorption/desorption  parameter  values. 
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Figure  5.43.  Comparison  of  the  tetrachloroethylene  plume  a)  mass  in 
solution,  and  b)  first  spatial  moment,  estimated  from  the 
data  and  simulated  by  an  equilibrium  transport  model  assum¬ 
ing  nonlinear  sorption/hysteretic  desorption  and  using 
experimentally  obtained  sorption/desorption  parameter 
values. 


126 


»  -  «  mM 


FIRST  MOMENT  (m)  ,N  SOLUTION  (fl) 


TIME  (doy») 


Figure  5.44.  Comparison  of  the  tetrachloroethylene  plume  a)  mass  in 
solution,  and  b)  first  spatial  moment,  estimated  from  the 
data  and  simulated  by  an  equilibrium  transport  model  assum¬ 
ing  nonlinear  sorption/hysteretic  desorption  and  using 
fitted  sorption/desorption  parameter  values* 
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Interestingly,  the  mathematical  form  of  this  model  is  identical  to  the 
physical  nonequilibrium  first-order  rate  model  (Nkedi-Kizza  et  al., 
1984).  Thus,  distinguishing  between  the  two  models  based  on  their  simu¬ 
lated  responses  is  impossible.  In  Chapter  4,  it  was  shown  that  an 
approximate  conversion  of  the  first-order  rate  model  to  the  diffusion 
models  is  possible,  and  that  attempting  to  distinguish  between  the  two 
types  of  models  is  difficult.  On  mathematical  grounds,  therefore,  chem¬ 
ical  nonequilibrium  is  as  adequate  an  explanation  of  the  Borden  observa¬ 
tions  as  is  physical  nonequilibrium. 

However,  if  we  consider  the  sorption  of  the  organic  solutes  in 
terms  of  the  widely  held  view  of  hydrophobic  partitioning  (Chiou  et  al., 
1979;  Karickhoff,  1984;  Gschwend  and  Wu,  1985),  it  would  seem  that  no 
chemical  reaction,  in  the  traditional  sense,  is  occurring,  and  that  slow 
sorption  is  really  slow  diffusion  to  a  sorption  "site"  rather  than  a 
slow  reaction  at  a  sorption  site. 

Biotransformation 

Roberts  et  al.  (1986)  noted  that  to  reproduce  the  observed  deceler¬ 
ation  with  time,  biotransformation  would  have  to  occur  selectively,  such 
that  solute  concentrations  at  the  leading  edge  of  the  plume  would  dimin¬ 
ish  more  rapidly  than  solute  concentrations  at  the  trailing  edge.  Current 
hypotheses  regarding  solute  transformation  in  groundwater  suggest  this 
pattern  is  unlikely,  and  in  fact  that  the  reverse  is  true:  biotransfor¬ 
mation  should  be  greater  at  the  trailing  edges  because  the  native  micro¬ 
organisms  capable  of  transforming  a  particular  compound  have  had  more 
time  to  multiply  or  acclimate  to  the  compound  (Bouwer  and  McCarty, 
1984).  Also,  the  organic  compounds  show  similar  decreases  in  velocity, 
regardless  of  whether  the  loss  in  mass  is  substantial  (bromoform)  or  not 
(carbon  tetrachloride).  For  these  reasons,  biotransformation  does  not 
seem  to  explain  adequately  the  observed  spatial  moment  behavior. 


External  Mass  Transfer 


The  physical  nonequilibrium  models  discussed  here  make  the  assump¬ 
tion  that  external  mass  transfer  resistance  is  negligible.  In  this  sec¬ 
tion,  the  validity  of  that  assumption  will  be  assessed. 

External  mass  transfer  resistance,  also  called  film  transfer  resis¬ 
tance,  is  due  to  a  boundary  layer  of  stagnant  water  surrounding  the 
solid  aquifer  particles.  For  sorption  to  occur,  solute  molecules  must 
first  diffuse  through  this  boundary  layer.  If  the  time  required  for  a 
solute  molecule  to  diffuse  through  this  layer  is  comparable  to  the  time 
required  for  the  molecule  to  diffuse  within  the  sorbing  particle  itself 
(internal  diffusion),  external  resistance  is  important  and  should  be 
accounted  for.  External  mass  transfer  is  typically  described  in  terms 
of  a  first-order  rate  expression,  with  kj  defined  as  the  first-order 
external  mass  transfer  coefficient. 

Crittenden  et  al .  (1986)  postulated  criteria  by  which  the  relative 
effects  of  dispersion,  internal  diffusion,  and  external  mass  transfer 
could  be  assessed.  According  to  Crittenden  et  al.  (1986),  when 
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the  influence  of  dispersion,  internal  diffusion,  and  external  mass 
transfer  on  spreading  will  be  equal.  If  one  or  two  of  the  three  terms 
in  Eq.  5-7  are  significantly  greater  than  the  other  terms,  the  mechanism 
represented  by  the  larger  term  or  terms  can  be  ignored. 

Table  3.8a  defines  an  effective  dispersion  coefficient  for  the 
diffusion  models  as: 


D  «  °x_  +  y26 
eff  1+6  a(l+B)3 


(5-8) 


The  equality  represented  by  the  first  two  terms  in  Eq.  5-7  is  the  same 
as  the'  equality  obtained  by  setting  the  two  terms  on  the  right-hand  side 
of  Eq.  5-8  equal  to  each  other.  That  is,  the  contributions  of  disper¬ 
sion  and  internal  diffusion  to  an  effective  dispersion  coefficient  are 
set  equal  to  each  other  by  the  criteria  of  Eq.  5-7.  Similarly,  the  last 


term  in  Eq.  5-7  represents  the  contribution  to  spreading  due  to  external 
mass  transfer  resistance.  Since  external  mass  transfer  is  described  by 
a  first-order  rate  expression,  a  first-order  rate  constant  a  may  be 
defined  as: 
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_ tn  t 
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Substituting  this  expression  into  Eq.  5-7  gives: 
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In  Table  3.8a  an  effective  dispersion  coefficient  for  the  first-order 
rate  model  is  defined  as: 
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Setting  the  two  terms  on  the  right-hand  side  of  Eq.  5-11  equal  to  each 
other  gives  the  equality  shown  between  the  first  and  third  terms  in 
Eq.  5-10.  Thus,  the  spreading  criteria  suggested  by  Crittenden  et  al. 
(1986)  based  on  inspection  of  the  transport  equation  solutions  for  uni¬ 
form,  one-dimensional  flow  are  seen  to  be  tantamount  to  equating  the 
second  moments.  Examination  of  Table  3.8b  reveals  that  Crittenden  et 
al.'s  (1986)  criteria,  derived  for  one-dimensional  flow,  are  also  applic¬ 
able  for  three-dimensional  flow  along  the  path  of  advective  transport. 

Using  Eq.  5-7  for  a  given  set  of  parameter  values  allows  the  deter¬ 
mination  of  the  relative  impact  of  dispersion,  internal  diffusion,  and 
external  mass  transfer  on  spreading.  Using  the  following  parameter 
values  for  tetrachloroethylene ,  the  compound  which  seems  to  be  most 
affected  by  physical  nonequilibrium  mechanisms: 


p  ■  1.81  g/cm3 
8  -  0.33 

Dmx  “  °'02  “2/d 
Kd  -  0.48  cm3/g 

d>  -  0.979 
vffl  *  0.091  ra/d 


d;  -  6.83  x  10'7  m2/d 
b  ■  0.001  m 
l  ■  0.29  m 
B  -  2.71 
f  -  0.0 
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gives  the  following  values  for  the  first  two  terms  of  Eq.  5-7: 
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D1- 

x 

*t(l  *  B>2  -  1.31 
Bv 

To  evaluate  the  third  term  in  Eq.  5-7,  an  estimate  for  kf  is  needed. 
Wilson  and  Geankopolis  (1966)  developed  the  following  correlation  to 
determine  kj  in  terms  of  the  Reynolds  (Re)  and  Schmidt  (Sc)  numbers: 


where: 
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Using  the  above  correlation,  and  the  value  of  DQ  for  tetrachloroethylene 
obtained  using  the  Wilke  and  Chang  (1955)  correlation,  gives: 


i 


a 


Sc  -  1265 

Re  -  2.11  x  10-3 

k^  •  0.115  m/d 

Using  this  value  of  kj  in  Eq.  5-7,  the  third  term  in  the  equation  is 
found  to  be: 
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From  a  comparison  of  the  three  terms  in  Eq.  5-7,  it  is  apparent 
that  although  the  impacts  on  spreading  due  to  dispersion  and  Internal 
diffusion  are  virtually  equal,  the  effect  due  to  external  mass  transfer 
resistance  is  three  orders  of  magnitude  less.  Thus,  it  seems  that  for 
any  reasonable  value  of  kf,  the  effect  of  external  mass  transfer  would 
be  negligible  compared  to  dispersion  and  internal  diffusion  effects. 
This  fact  is  illustrated  in  Figure  5.45.  In  Figure  5.45,  an  orthogonal 
collocation  code  developed  by  Crittenden  et  al.  (1986),  which  combines 
one-dimensional  advective/dlspersive  transport  with  Internal  diffusion 
into  spheres  and  external  mass  transfer,  is  used  to  obtain  breakthrough 
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Figure  5.45.  Simulations  of  breakthrough  responses  for  various  values 
of  the  external  film  transfer  coefficient. 


responses  using  the  parameter  values  for  tetrachloroethylene .  The  fig¬ 
ure  demonstrates  that  varying  kf  over  four  orders  of  magnitude  (i.e., 
increasing  or  decreasing  by  two  orders  of  magnitude  relative  to  the  best 
estimate  of  kf  *  0.115  m/d)  has  no  effect  on  the  response.  Thus,  it  is 
clear  that  the  effect  of  external  mass  transfer  resistance  is  negligi¬ 
ble,  for  the  parameter  values  under  consideration. 


Aquifer  Heterogeneity 

Roberts  et  al.  (1986)  postulated  that  a  trend  of  increasing  sorp¬ 
tion  in  the  direction  of  plume  movement  might  lead  to  the  observed  loss 
of  solute  mass  in  solution  and  the  deceleration  of  the  solute  plumes. 

Durant  (1986)  measured  the  uptake  of  tetrachloroethylene  by  Borden 
aquifer  material.  The  aquifer  material  was  obtained  from  cores  taken 
along  the  path  of  plume  movement*  A  linear  best  fit  of  the  measured 
values  of  tetrachloroethylene 's  distribution  coefficient  (K^)  versus 
plume  travel  distance  along  the  path  of  tetrachloroethylene  plume  move¬ 
ment  (i*e.,  cores  S2-S5  of  Durant,  1986)  indicates  an  initial  K.,  at  the 


area  of  solute  Injection,  of  0.431  cor/g,  with  an  Increase  of  0.016  cnr/g 
per  meter  of  plume  travel  distance.  Assuming  local  equilibrium,  these 
values  of  Kj  can  be  converted  to  corresponding  values  of  a  retardation 
factor  (R).  Performing  this  conversion  gives  an  Increase  in  R  from 
about  3.4  at  the  area  of  plume  injection,  to  about  4.5  at  a  point  13 
meters  dovngradient,  where  the  center  of  mass  of  the  tetrachloroethylene 
plume  was  located  at  the  last  synoptic  sampling  session.  This  increase 
in  R  is  significantly  less  than  the  Increase  in  R  from  3  to  6  shown  in 
Figure  8  of  Roberts  et  al.  (1986).  The  increase  of  3  to  6  was  estimated 
based  on  the  deceleration  of  the  tetrachloroethylene  plume  over  the 
course  of  the  experiment.  A  similarly  large  increase  in  R  would  be  esti¬ 
mated  based  on  the  decrease  in  the  total  mass  of  tetrachloroethylene  in 
solution  during  the  experiment. 

Thus,  although  a  trend  of  increasing  sorption  in  the  direction  of 
plume  movement  may  contribute  to  the  observed  loss  in  mass  and  decelera¬ 
tion  of  the  solute  plumes,  the  extent  of  this  contribution  appears  to  be 
small  when  compared  with  the  magnitude  of  the  observed  effect. 

CONCLUSIONS 

It  has  been  shown  that  concentration  distribution  histories  ob¬ 
served  during  a  natural  gradient  experiment  studying  solute  transport 
may  be  Interpreted  using  a  physical  nonequilibrium  model  which  hy¬ 
pothesizes  advective/di8persive  solute  transport  combined  with  solute 
diffusion  into  regions  of  immobile  water.  Parameter  values  obtained 
independently  were  shown  inadequate  to  model  spatial  moment  behavior. 
However,  qualitative  aspects  of  the  spatial  data  were  simulated  by 
assuming  hindered  diffusion  within  spherical  Immobile  regions,  or  a  high 
sorption  capacity  within  layered  regions.  The  validity  of  these  assump¬ 
tions  is  currently  under  experimental  investigation.  Making  the  assump¬ 
tion  of  diffusion  into  layers  of  high  sorption  capacity,  the  model  suc¬ 
cessfully  simulated  temporal  behavior,  using  parameters  obtained  from 
spatial  data  only,  thus  providing  an  example  of  the  complementarity  of 
spatial  and  temporal  behavior.  This  complementarity  was  also  demon¬ 
strated  by  comparing  results  of  spatial  and  temporal  moment  analyses  at 
comparable  time/distance  scales.  General  transport  behavior  inferred 
from  the  spatial  data  could  be  inferred  from  the  temporal  data  also. 


Although  a  model  based  on  the  traditional  assumptions  of  equilib¬ 
rium  advective/dispersive  transport  adequately  describes  many  aspects  of 
the  observed  behavior,  several  phenomena,  such  as  the  decline  in  organic 
solute  mass  in  solution,  the  deceleration  of  the  sorbing  solute  plumes, 
and  the  tailing  of  the  near-field  breakthrough  responses,  are  not 
consistent  with  the  equilibrium  model. 

Other  models  were  also  examined.  In  particular,  numerical  simula¬ 
tions  were  performed  to  assess  the  spatial  moment  behavior  of  solute 
concentration  distributions  predicted  by  models  which  assume  nonlinear 
and/or  hysteretic  sorption.  It  was  discovered  that  nonlinear  sorption 
could  also  qualitatively  explain  the  spatial  and  temporal  moment  behav¬ 
ior,  though  at  the  low  concentrations  under  consideration,  there  was 
little  experimental  evidence  of  nonlinear  sorption  behavior.  Neverthe¬ 
less,  nonlinear  sorption  does  offer  a  possible  explanation  for  the 
observed  transport  behavior  as  well. 
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CHAPTER  6 

CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 

The  focus  of  this  work  has  been  the  development,  analysis,  and 
application  of  a  three-dimensional  solute  transport  model  which  incorpo¬ 
rates  the  mechanism  of  diffusion  into  the  traditional  advection/disper- 
sion  transport  equation.  In  this  chapter,  conclusions  drawn  from  this 
study  will  be  presented,  and  recommendations  made  for  related  future 
work. 

The  specific  conclusions  from  this  study  can  be  summarized  as  fol¬ 
lows: 

1.  Analytical  solutions  to  various  formulations  of  three-dimensional 
physical  nonequilibrium  solute  transport  models  were  presented. 
The  solutions  approach  the  well-known  solution  of  the  three- 
dimensional  advectlon/dispersion  equation  In  the  limiting  case 
where  physical  nonequilibrium  is  negligible. 

2.  Expressions  for  the  spatial  and  temporal  moments  of  concentration 
distributions  simulated  by  the  three-dimensional  physical  nonequi¬ 
librium  models  were  derived,  using  a  modified  form  of  Aris'  method 
of  moments.  Analysis  of  the  moments  showed: 

a.  The  zeroth  spatial  moment  of  the  concentration  dis¬ 
tribution  of  solute  in  the  mobile  water  region  is  a  decreasing 
function  of  sampling  time.  Thus,  the  physical  nonequilibrium 
models  predict  a  loss  of  solute  mass  In  the  mobile  region 
over  time. 

b.  Although  the  zeroth  and  first  temporal  moments  are 
independent  of  the  rate  of  solute  transfer  Into  the  immobile 
water  region,  all  the  spatial  moments  are  dependent  on  the 
mass  transfer  rate.  One  implication  of  this  is  that  while 
the  retardation  factor  obtained  from  temporal  breakthrough 
response  data  is  Independent  of  sampling  location,  the  retar¬ 
dation  factor  obtained  from  synoptic  data  is  an  Increasing 
function  of  sampling  time. 

c.  In  previous  research,  expressions  for  equilibrium 
model  parameters  which  give  responses  equivalent  to  physical 
nonequilibrium  models  have  been  derived  using  temporal  moment 
analyses.  This  work  shows  that  the  long  time  values  of  the 


equivalent  parameters  obtained  using  spatial  moment  analyses 
are  the  same  as  the  values  obtained  using  temporal  moments. 

Thus,  the  expressions  derived  previously,  based  on  temporal 
distributions,  are  also  applicable  when  considering  spatial 
distributions. 

3.  Due  to  the  effect  of  one-dimensional  versus  three-dimensional  dis¬ 
persion,  the  first  temporal  moment  obtained  using  a  one-dimensional 
model  is  greater  than  the  first  temporal  moment  obtained  using  a 
three-dimensional  model  given  the  same  velocity.  This  is  true  for 
both  equilibrium  and  nonequilibrium  models. 

4.  At  high  values  of  Pem  (Pem  >  50),  the  effect  of  the  boundary  condi¬ 
tions  on  temporal  simulations  of  the  models  becomes  small. 

5.  Based  on  the  equality  of  the  first  and  second  moments,  equivalent 
parameters  were  defined  for  equilibrium,  first-order  rate,  and 
diffusion  models.  The  differences  between  the  model  simulations 
using  these  equivalent  parameters  are  not  large,  and  differentiat¬ 
ing  among  these  models  based  on  measured  breakthrough  response  data 
would  be  difficult,  owing  to  parameter  and  measurement  uncertainty. 
Differentiating  between  equilibrium  and  nonequilibrium  models  may 
be  possible  using  spatial  moment  data,  by  taking  advantage  of 
qualitative  differences  between  the  spatial  moment  predictions  of 
the  two  types  of  models. 

6.  Equilibrium  and  nonequilibrium  models  were  applied  to  data  obtained 
in  a  field  experiment  studying  solute  transport  in  an  unconfined, 
sand  aquifer.  Although  the  experiment  was  conducted  in  a  rela¬ 
tively  homogeneous  aquifer,  and  indeed,  the  equilibrium  model 
adequately  simulated  much  of  what  was  observed,  some  aspects  of  the 
solute  behavior  could  not  be  explained  using  the  traditional  equi¬ 
librium  approach.  The  assumption  of  physical  nonequilibrium  offers 
a  reasonable,  simple  explanation  which  qualitatively,  and  to  some 
extent  quantitatively,  explains  experimental  observations. 

The  following  topics  are  recommended  for  future  research: 

1.  Development  of  methods  to  measure  independently  such  parameter  val¬ 
ues  as  the  diffusion  rate  constant  and  the  ratio  of  mobile  to  total 
water  content  is  crucial  both  for  gaining  further  understanding  of 
the  diffusion  mechanism  and  for  using  the  models  predlctively. 
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Further  improvement  in  the  numerical  method  used  to  evaluate  first 
and  second  spatial  moments  is  needed.  The  method  presented,  though 
useful,  is  cumbersome. 

The  analytical  solutions  presented  in  this  study  are  obviously 
limited  in  their  application  to  real  field  situations,  where  condi¬ 
tions  vary  both  spatially  and  temporally.  The  development  of  a 
three-dimensional  numerical  transport  model  which  incorporates 
physical  nonequilibrium  mechanisms  (using  the  analytical  solutions 
presented  here  as  developmental  tools)  would  allow  these  physical 
nonequilibrium  models  to  be  generally  applied  to  field  transport 
situations. 

The  field  experiment  discussed  in  this  study  was  conducted  in  a 
relatively  homogeneous,  sandy  aquifer,  where  physical  nonequi¬ 
librium  effects  were  not  expected  to  play  an  important  role. 
Nevertheless,  in  this  work,  considerable  evidence  is  presented 
indicating  that  physical  nonequilibrium  mechanisms  may  have 
impacted  the  solute  transport.  Conduct  of  a  field  experiment  in  an 
aquifer  where  physical  nonequilibrium  is  expected  to  play  a  more 
dominant  role  (e.g.,  an  aggregated  or  stratified  medium)  may  shed 
more  light  on  how  this  mechanism  affects  transport. 

Based  on  the  usefulness  of  the  spatial  data  obtained  from  the 
Borden  field  experiment,  it  seems  likely  that  further  field  studies 
will  be  conducted  which  attempt  to  produce  corroborative  data. 
The  analysis  of  other  transport  models  with  respect  to  three- 
dimensional  spatial  concentration  distributions  (similar  to  the 
limited  one-dimensional  analysis  conducted  in  Chapter  5  for  the 
hysteretic  and  nonlinear  sorption  models)  would  assist  in  analyzing 
spatial  data  obtained  from  the  field. 


APPENDIX  A 

DERIVATION  OF  THE  SOLUTION  TO  THE  FIRST-ORDER  RATE  MODEL 
FOR  AN  INSTANTANEOUS  POINT  SOURCE 


Differential  Equations 
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Initial/Boundary  Conditions 

Cm(x,y,z,0)  -  M36(x)6(y)6(z) 

(A-3a) 

Cm(±»,y,z,t)  -  Cm(x,±»,z,t)  -  Cm(x,y,±«,t)  -  0 

(A-3b) 

ctm(x,y,z,°)  -  0 

(A-3c) 

Basically  following  the  methods  of  Lindstrom  and  Narasimhan  (1973), 
who  solved  a  similar  set  of  equations,  take  the  Laplace  transform  of 

(A-2),  and  then  apply  (A-3c)  to  obtain: 
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where  the  Laplace  transform,  F(x,y,z,s),  of  a  function,  F(x,y,z,t),  Is 
defined  as: 
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F(x,y,z,s)  -  /  e  stF(x,y,z,t)dt  (A-5) 
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Now,  find  the  Laplace  transform  of  (A-l),  using  (A-3a)  and  (A-A)  to 
obtain  the  following  differential  equation: 
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with  the  boundary  condition: 


C  (±»,y,z,s)  -  C  (x,±*,z,s)  -  C  (x,y,±«,s)  *  0 
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(A-7) 


A  solution  to  (A-6)  and  (A-7)  nay  be  found  using  Fourier  transforms. 
Define  F(p,y,z,s)  as  the  Fourier  transform  of  F(x,y,z,s)  such  that: 


£(p,y.z,s)  ■  /  F(x,y,z,s)e  lpxdx 
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where  p  is  the  Fourier  transform  variable  in  the  x-direction.  Similarly, 
define  q  and  u  as  the  Fourier  transform  variables  in  the  y-  and 
z-directions,  respectively.  Taking  the  Fourier  transform  of  (A-6)  in 
the  x-,  y-,  and  z-directions,  yields: 
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To  obtain  (^(x.y.z.s),  it  is  necessary  to  invert  the  Fourier  transform. 
In  the  x-direction  apply  the  equation: 


F(x,y,z,s)  m  j  £(P»y»*»s>  eipXdp 


(A-10) 


Analogous  equations  will  be  used  to  invert  the  Fourier  transformed  solu¬ 
tion  in  the  y-  and  z-directions. 

Simplifying  (A-9)  to  read: 
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and  using  the  definite  integrals  (Gradshteyn  and  Ryzhik,  1980): 
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for  x  >  0,  it  is  found  that  the  Fourier  inverse  of  (A-9)  in  the  x-direction 
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To  obtain  the  inverse  Fourier  transform  of  (A-13)  in  the  y-direction, 
simplify  (A-13)  to  read: 
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Applying  the  y-direction  analog  of  Eq.  A-10  and  using  the  definite  inte¬ 
gral  (Gradshteyn  and  Ryzhik,  1980): 
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gives 
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To  obtain  the  inverse  Fourier  transform  of  (A-16)  in  the  z-direction, 
once  again  follow  the  same  basic  method.  Simplify  (A-16)  to  read 
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Apply  the  z-direction  analog  of  (A-10),  using  the  definite  Integral 
(Gradshteyn  and  Ryzhik,  1980) 
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The  inverse  Laplace  transform  of  (A-19)  can  be  obtained  using  the 


methods  of  Lindstrom  and  Narasimhan  (1973).  The  final  result  is 
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APPENDIX  B 

DERIVATION  OF  THE  SOLUTION  TO  THE  DIFFUSION  MODEL 
WITH  SPHERICAL  IMMOBILE  ZONES 


Differential  Equations 
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Initial/Boundary  Conditions 
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Following  the  methods  of  van  Genuchten  et  al.,  (1984),  take  the 
Laplace  transform  of  (B-2),  applying  Initial  condition  (B-4c)  to  obtain 
the  following  ordinary  differential  equation: 
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Transforming  boundary  conditions  (B-4d)  and  (B-4e)  gives 


Solve  the  above  ordinary  differential  equation  with  boundary  conditions 
to  find 
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Now  take  the  Laplace  transform  of  (B-3) 
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and  substitute  (B-7)  into  (B-8)  to  obtain 
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Take  the  Laplace  transform  of  (B-l),  using  (B-9)  and  Initial  conditions 
(B-4a)  and  (B-4c)  to  get  the  following  differential  equation 
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Transforming  boundary  condition  (B-4b)  gives 


Cm^i®,y,z*8^  *  Cm(x,±co,z,s)  -  Cm(x,y,dt«,s)  -  0 
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A  solution  to  (B-10)  and  (B-ll)  may  be  found  using  Fourier  transforms. 
As  in  Appendix  A,  the  Fourier  transformed  solution  is 
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Note  the  similarity  between  Eq.  B-12  and  Eq.  A-9.  To  obtain  Cm(x>y»z>s), 
the  methods  of  inverting  the  Fourier  transform  in  the  x-,  y-,  and  z- 
directions  used  in  Appendix  A  are  directly  applicable.  The  result  is 


APPENDIX  C 

DERIVATION  OF  ABSOLUTE  SPATIAL  MOMENTS 
FOR  THREE  SOLUTE  TRANSPORT  MODELS 


Local  Equilibrium  Model 
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Taking  the  Fourier  transform  of  the  equation  and  initial  condition, 
where  p,  q,  u  are  the  Fourier  transform  variables  in  the  x-,  y-, 

z-directions,  respectively,  yields: 
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Applying  the  expression 
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requires  the  differentiation  of  C(p,q,a,t)  with  respect  to  p,  q,  and  u. 
This  operation  is  easily  done*  Taking  the  limit  as  the  transform  vari¬ 
ables  approach  0,  gives: 
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First-Order  Rate  Model 
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Following  Appendix  A,  these  equations  may  be  solved  in  the  Laplace  and 
Fourier  domains,  to  obtain 
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(see  Eq.  A-9,  Appendix  A). 

The  inverse  Laplace  transform  of  this  equation  needs  to  be  derived. 
Simplify,  by  defining 
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Multiply  the  numerator  and  denominator  of  the  right-hand  side  of  Gq. 
C-12  by  8  +  a,  to  obtain 
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Defining  G  and  F  as  the  roots  of  the  quadratic  equation  in  the  denomina¬ 
tor  of  the  right-hand  side  of  Gq.  C-13  gives: 
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where  E  and  F  are  functions  of  A,  a,  and  g. 

Differentiation  of  Gq.  C-14  is  straightforward,  though  tedious 
Apply  Gq.  C-6,  to  find: 
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Diffusion  Models 
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(C-15) 


(C-16) 

(C-17) 

(C-18) 


(C-19a) 

(C-19b) 


(^(i-.y.z.t)  -  (^(x.i-.z.t)  -  (^(x.y.i-.t)  -  0 
Ca(0,x,y,z,t)  +  • 

C8(b,x,y,z,t)  -  (^(x.y.z.t) 


(C-19c) 


(C-19d) 

(C-19e) 


Following  Appendix  B,  these  equations  may  be  solved  in  the  Laplace 
and  Fourier  domains,  to  obtain 


M, 


Cm(p,q,u,s)  -  - j 


where 


Dxp“  +  Dyq2  +  Dzu2  +  vip  +  N2 


36s  i^(uib) 

+  s  for  the  spherical  diffusion  model 


(C-20) 


2fis  I..(iijb) 
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ub  I  (ub) 
o 


s  for  the  cylindrical  diffusion  model 


!ln!l  ‘‘’c  +  8  for  the  layered  diffusion  model 
ub  Cosh  tub 


To  derive  the  inverse  Laplace  transform  of  Eq.  C-20,  it  is  neces¬ 
sary  to  apply  Bromwich's  complex  inversion  formula: 

.  a+i» 

™  2ui  /  e  ^m 
a-i® 

Using  the  same  technique  as  Rosen  (1954),  Rasmuson  and  Neretnieks  (1980), 
and  van  Genuchten  et  al.  (1984),  who  derived  similar  inversions,  it  can 
be  shown  that 


«  M«  -ie  i®  ie  st 

Cm(p,q,u,t)  -  lim  {  /  +  /  +  /  )  fiTSy  d8  (C"21) 

e+0  -i®  ie  -ie 

where  n(s)  ■  DXP2  +  Dyq^  +  Dzu*  +  vip  +  N^. 

Perform  the  above  integration,  basically  following  the  methods  of 
Rasmuson  and  Neretnieks  (1980)  and  van  Genuchten  et  al.  (1984),  to  find: 
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b 
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D  D 

and  ^  and  ^  are  as  defined  in  Chapter  2  for  the  different  immobile 
region  geometries.  As  p,  q,  and  u  approach  zero,  and  ate  as 
defined  in  Appendix  B  (with  v  *  0). 

Use  Eq.  C-6,  to  write 
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(C-22) 


Differentiating  Eq.  C-22  by  applying  Liebnitz’  rule,  and  then  applying 
Eq.  C-6,  leads  to  the  following  expressions: 
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where  for  spherical  diffusion: 
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for  cylindrical  diffusion: 


for  layered  diffusion: 
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APPENDIX  D 

DERIVATION  OF  ABSOLUTE  SPATIAL  MOMENTS  OF  THE 
IMMOBILE  REGION  SOLUTE  DISTRIBUTION 


First-Order  Rate  Model 


In  Appendix  A  it  was  found 


Clm(x,y,z,s)  -  Cm(x,y,Z,s) 


(D-l) 


Take  the  Fourier  transform  of  (D-l): 
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In  Appendix  C  it  was  found 
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Combining  (D-2)  and  (D-3): 
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Clni(p,q,u,8)  -  [(8  +  0)  -  E][(s  +  a)  -  F] 
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Follow  Appendix  C  to  invert  the  Laplace  transform,  to  obtain 
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Finally,  differentiate  Eq.  D-5  and  apply  Eq.  3-23  to  derive  the  follow¬ 
ing  absolute  moments: 
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Follow  Appendix  B,  to  find 


C.  (x.y.z.s)  -  8*  C  (x,y,z,s) 


sB  m 


(D-6) 


where 


3Bs  i,((ob) 

Nz  -  — ,  rr  + 


uib  i  (o)b) 
o 


for  the  spherical  diffusion  model 


28s  I .(mb) 
-  - - - - 


u>b  I  (t»)b)  8 

o 


for  the  cylindrical  diffusion  model 
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Take  the  Fourier  transform  of  (D-6): 


Cim(p,q,u,s)  -  Cm(p,q,u,s) 


(D-7) 


In  Appendix  B  it  was  found 
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Combining  (D-7)  and  (D-8)  yields: 
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Defining: 
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gives: 
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Follow  Appendix  C  to  obtain  the  inverse  Laplace  transform  of  Cjm(p,q,u,s) , 
using  Bromwich's  complex  inversion  formula,  to  find: 
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where  A  -  Dxp2  +  Dyq2  +  Dzu2  +  vlp 


v( v  +  2)D 


v  ■  3  for  spherical  diffusion 
where  v  ■  2  for  cylindrical  diffusion 
v  *  1  for  layered  diffusion 


and  and  122  are  defined  in  Appendix  C.  Differentiating  and  applying 
Eq.  3-23  leads  to: 
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APPENDIX  E 

BREAKTHROUGH  RESPONSE  DATA 
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Figure  E.3.  Breakthrough  response  at  well  location  x  *  5.00  m, 
y  *  0.00  m,  and  z  -  -3.26  m. 


Figure  E.4.  Tetrachloroethylene  breakthrough  response  at  well 
location  x  ■  10.00  m,  y  *  4.60  m,  and  z  *  -3.88  m. 
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Figure  E.7.  Tetrachloroethylene  breakthrough  response  at  well 
location  x  -  13.10  m,  y  -  4.05  m,  and  z  *  -3.72  m. 


NORMALIZED  CONCENTRATION  ►*  NORMALIZED  CONCENTRATION 


jure  E.ll.  Breakthrough  response  at  well  location  x  -  21.00  m 
y  »  9.00  m,  and  z  *  -4.77  m. 


Figure  E.12.  Bromide  breakthrough  response  at  well  location  x 
24.00  o,  y  *  9.00  m,  and  z  ■  -4.76  m. 


APPENDIX  F 

HYSTERESIS  EXPERIMENT  METHODOLOGY 


A  bottle- pint  partitioning  technique  to  measure  long-term  sorption 
(Ball  and  Roberts,  1985)  was  combined  with  a  technique  to  measure  long¬ 
term  desorption.  In  order  to  assess  whether  sorption/desorption  of 
tetrachloroethylene  onto  Borden  aquifer  material  was  hysteretic. 

Sorption  Experiments 

Sorption  experiments  were  conducted  using  the  bottle-point  parti¬ 
tioning  technique  of  Ball  and  Roberts  (1985)  to  measure  long-term  sorp¬ 
tion  while  carefully  controlling  losses.  Using  this  technique,  a  known 
mass  of  groundwater  synthesized  to  approximate  aquifer  conditions 
(Uurtis  et  al.,  1986)  and  a  known  mass  of  sorbent,  consisting  of  homoge¬ 
nized  Borden  aquifer  material,  were  added  to  glass  ampules  and  the 
sorbent  allowed  to  come  to  complete  water  saturation.  A  known  amount  of 
carbon-14  labeled  solute  was  then  Injected  into  the  glass  ampule,  and 
the  ampule  Immediately  flame-sealed.  After  sealing,  the  sample  was 
allowed  to  equilibrate  over  a  30-day  period.  After  the  30  days,  half 
the  samples  were  centrifuged  (2000  rpm  for  30  min),  cracked  open,  and 
supernatant  immediately  transferred  to  scintillation  cocktail  for  count¬ 
ing  on  a  Packard  Tricarb  Model  4530  Scintillation  Counter.  This  count 
provided  an  accurate  measure  of  the  solute  concentration  in  the  aqueous 
phase.  By  difference,  the  sorbed- fhase  concentration  was  determined, 
thus  determining  a  point  on  the  sorption  isotherm.  Losses  due  to  parti¬ 
tioning  to  the  glass  ampule  and  to  the  small  volume  of  air  in  the  ampule 
were  accurately  and  reproducibly  quantified  using  blanks  without  soil. 

Desorption  Experiments 

Desorption  experiments  were  then  conducted  upon  those  samples  not 
opened  in  the  sorption  experiments.  These  samples  were  placed  into 
specially  fabricated  stainless  steel  vessels.  To  minimize  losses,  the 
vessels  were  constructed  to  be  airtight  and  to  have  only  stainless  steel 
in  contact  with  the  sample.  Synthetic  groundwater,  which  had  been  pre- 
equillbrated  with  the  aquifer  material,  was  added  to  the  vessel,  the  cap 
to  the  vessel  closed,  and  the  vessel  shaken  vigorously  and  centrifuged 
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for  30  min.  The  shaking  and  centrifugation  broke  the  glass  ampule  in 
the  steel  vessel,  causing  the  aquifer  material  and  solute  within  the 
ampule  to  be  diluted  by  the  synthetic  groundwater  in  the  vessel.  The 
sample  was  then  allowed  to  equilibrate  for  30  days.  After  the  30  days, 
the  steel  vessel  was  centrifuged,  opened,  and  the  supernatant  trans¬ 
ferred  to  scintillation  cocktail  for  counting.  As  in  the  sorption 
experiment,  aqueous-phase  concentration  was  thus  directly  measured,  and 
the  sorbed-phase  concentration  determined  by  difference.  Blanks  showed 
that  losses  over  the  course  of  the  60-day  sorption/desorption  experiment 
amounted  to  less  than  15%. 
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